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Final  Report 

Turbulence  Modeling  using  Body  Force  Potentials. 

Objective 

The  objective  of  this  study  was  the  development  and  testing  of  the  turbulent  potential  model  in  complex 
non-equilibrium  flows.  The  initial  emphasis  was  on  evaluating  and  enhancing  the  models  performance  in 
rotating  flows  and  in  transition  to  turbulence.  The  project  was  extended  in  the  third  year  to  examine  the 
model’s  ability  to  act  as  a  subgrid  scale  model  in  Large  Eddy  Simulations  (LES). 

The  turbulent  potential  model  is  an  alternative  approach  to  Reynolds  Averaged  Navier-Stokes  (RANS) 
turbulence  modeling,  where  the  primary  quantity  of  interest  is  no  longer  the  Reynolds  stress  tensor.  Instead 
the  divergence  of  the  Reynolds  stress  tensor,  a  body  force  vector,  is  the  primary  quantity  describing  how 
the  turbulence  affects  the  mean  flow  evolution,  and  transport  equations  for  the  Helmholtz  potentials  of  that 
body  force  vector  are  derived  and  then  closed  using  modeled  transport  equations. 

The  potential  model  approach  does  not  hypothesize  any  explicit  relationship  between  the  turbulence  and 
the  mean  flow,  it  is  therefore  capable  of  capturing  non-equilibrium  turbulent  flows,  and  is  physically 
equivalent  to  Reynolds  stress  transport  (RST)  models.  However,  because  the  focus  is  on  a  vector  quantity 
rather  than  a  tensor,  the  model  equations  for  the  turbulent  potential  model  are  simpler  than  RST  models  and 
roughly  comparable  in  cost  and  complexity  to  the  widely  used  two  equation  models. 


1.  Potential  Model 

1.1  Turbulence  Body  Force  Potentials 

Turbulence  potentials  were  introduced  by  Perot  (1996),  and  are  related  to  the  Reynolds  stress  tensor,  R, 
through  the  following  relationship 

V  R  =  Vx\|/  +  V(j)  (1.1) 

\J/  is  a  vector,  and  (f>  is  a  scalar.  They  can  be  regarded  as  potentials  because  their  curl  and  gradient  give  a 
body  force.  The  following  relationship  is  chosen  to  ensure  that  the  potentials  are  unique 

V-\|f  =  0  (2.1) 

The  physical  meaning  of  the  turbulent  potentials  is  explained  by  Perot  and  Moin  (1996).  <))  can  be 
interpreted  as  a  pressure  created  due  to  the  turbulence.  The  vector  potential  can  be  interpreted  as  a  the 
average  vorticity  in  turbulent  filaments.  If  the  eddy  viscosity  hypothesis  is  used  then  \| f  can  be  shown  to 
be  equal  to  an  eddy  viscosity  multiplied  by  the  mean  vorticity.  In  the  actual  model,  the  eddy  viscosity 
hypothesis  is  not  used  and  a  transport  equation  for  the  vector  potential  is  solved. 
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The  turbulent  potentials  are  related  to  R  by  the  following  equations 


V2(()  =  V-(V-R)  (1.3) 

V2\jf  =  -V  x  (V  •  R)  (1.4) 

which  can  be  obtained  by  operating  on  (1.1)  with  a  divergence  and  with  a  curl  respectively.  However,  if  the 
flow  is  inhomogeneous  in  only  one  direction  in  a  Cartesian  coordinate  system,  then  we  can  obtain  a  one-to- 
one  relationship  between  the  relevant  components  of  the  turbulent  potentials  and  the  Reynolds  stresses. 
Thus  if  the  flow  is  inhomogeneous  in  only  the  y  direction  (x2  in  Cartesian  tensor  notation),  then  (1.3)  and 
(1.4)  will  reduce  to 


^,22  —  R-22,22 

Vl,22  =  ~ ^23,22 

¥»=«  <l7> 

¥3, 22  =  *12,22 


Now,  if  we  assume  that  the  turbulent  potentials  vanish  when  the  turbulence  vanishes,  then  we  have,  from 
(1.5)  to  (1.8), 


<t>=R22  >  ¥,=-*32  J¥2=°J¥3  =  R,2  0-9> 

Flows  with  inhomogeneity  in  one  direction  include  fully  developed  channel  flows,  fully  developed 
spanwise  rotating  channel  flows,  fully  developed  streamwise  rotating  channel  flows,  and  (to  a  some  degree) 
2D  and  3D  boundary  layers  and  free  shear  flows. 


1.2  The  Potential  Model  Transport  Equations 

The  potential  model  is  a  set  of  transport  equations  for  the  turbulent  potentials  and  two  auxiliary  variables  (k 
and  e)  used  to  model  the  source  terms  in  those  transport  equations.  The  auxiliary  variables  are  not  used  to 
model  the  turbulence  effects  on  the  mean  flow  so  this  model  is  not  close  to  a  standard  k/e  model.  The 
equations,  their  derivation,  and  the  rational  for  the  model  is  an  extensive  topic,  that  was  completed  prior  to 
this  project  and  presented  in  Perot  (1999).  This  work  also  uses  the  model  to  calculate  rotating  channel  flow, 
a  mixing  layer,  adverse  pressure  gradient  boundary  layer  and  the  flow  over  a  backward  facing  step.  The 
results  are  all  in  good  accord  with  corresponding  DNS,  LES  or  experimental  data.  Some  more  practical 
engineering  flows  like  the  flow  due  to  an  impinging  jet  were  successfully  simulated  by  Wang  (2000).  The 
applications  of  the  model  to  3D  boundary  layer  and  rotating  turbulent  flows  were  studied  by  Are  (2001) 
and  Bhattacharya  (2002)  respectively.  More  recently,  Zhang  (2002)  applied  the  model  with  an  unstructured 
staggered  mesh  scheme  to  simulate  the  unsteady  vortex  shedding  in  flows  over  two-dimensional  obstacles 
and  Wang  (2002)  has  evaluated  the  models  ability  to  predict  transition.  The  work  of  Wang  and 
Bhattacharya  were  funded  by  this  project  and  is  presented  below. 

The  turbulent  potential  model  is  given  by  the  following  equations. 
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Dtp 

Dt 


■■v-{v+vt)v<p  +utw+nr;pld-e,+c^\-L5^y 


kVJ 


^ = v  ■  {v+v^vy+n^+n;^  -e¥ + c, 


[±Xcok+k-A 


(1.10) 


(1.11) 


Dk 

~  =  V-(v  +  v,ak)Vk  +  P-e 


(1.12) 


^  =  V  •  (v+ vla£)Ve+ Cr  UcE,P-Cne) 

Dt  k 

Constants  and  parameters  are  given  by: 

p=r-t.  v,=o.2i^-,  a=-\j,  ^=h4,  q=-2r(2a-m-onr) 

S  1+2  k  £ 


(1.13) 


£  =  £ /[1  +  lOv^l  Vk  |  Ik] 


V, 


C„=2.0— ^T’^00033’  = 


v,  +10t> 
Pa 


(1  +  0.105£) 
(l  +  0.06£> 


C;2  =  0.6  +  0.3—  + 1 .4(2#  - 1 .)  C;2  =  f  -  (0.6  -  C^2 )  +  C^2 


(1.14) 


ak=0.33+0.67P/i 
ct£=0.33  +  0.5P/£ 


( 


C  .  =  1 .45,  C  2  =1.83-0.16  exp 


„2\ 


-0.25- 


V£ 


V  ^  J 

A  =  0.02s  x  (s  x  \y)  +  0. 1 5(d)-  s)  x  ((d)  -  s)  x  \y) 

The  slow  pressure  strain  models  are 


nr  =  -(v+2t>,)V|  f  [Vk  +  C'.jila-W 


+  {Cp2+Cp^)- 


.-P/. 


v, 


<!> 


(1.15) 
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(1.16) 


nr 


=  -{v+2v,)V\¥-j-Vk-Cp]j(l-a)ys 


p]k 


The  fast  pressure  strain  models  are 

117"  =  (2  -  csp2  )(P  -  J<P) - (2  -  Cwp2  )(W  ■  d>) 
n;pid  ={2-C%){\-q)(t>si -(2 -C?3)q(M)-U 

The  dissipation  terms  are  given  by 


=  2vV<f>l/2  •  V^2  +  2a—  e 
k 


(1.17) 


The  vector  A  appears  only  in  3-D  flows,  and  therefore  contributes  in  predicting  the  cases  of  the  swirling 
pipe  flow  and  streamwise  rotating  channel  flow.  The  terms  involving  C,  are  the  transition  terms.  With 
typical  values  of  this  constant  these  terms  have  no  affect  on  fully  turbulent  flows.  In  the  presentation 
above,  the  rapid  pressure  terms  also  contain  the  production  terms  since  they  are  closely  related. 


2.  Rotation 
2.1  Background 

Rotation  occurs  in  many  turbulent  flows,  and  is  a  critical  aspect  of  such  devices  and  flows  as  pumps, 
turbines,  and  trailing  vortices.  Turbulence  is  insensitive  to  linear  accelerations  but  not  rotational  ones  (a 
spinning  top  has  the  same  transformation  characteristics).  Centrifugal  and  Coriolis  effects  have  both  direct 
and  indirect  effects  on  the  turbulence  evolution.  The  complexity  of  how  rotation  influences  turbulence 
evolution  can  be  seen  even  in  the  simplest  turbulent  flows.  The  decay  of  isotropic  turbulence  is  inhibited 
when  the  fluid  is  strongly  rotated.  Consequently,  high  turbulence  levels  can  persist  for  a  much  longer  time 
when  the  fluid  is  spinning  (such  as  in  the  core  of  a  trailing  vortex).  However,  while  the  turbulent  kinetic 
energy  remains  large,  the  influence  of  the  turbulence  on  the  mean  flow  (shear  stress)  actually  decreases 
compared  to  a  non-rotating  flow.  The  core  of  trailing  vortex  behaves  as  if  it  was  laminar  (not  influenced  by 
turbulence). 

These  experimental  results  indicate  that  simple  turbulence  models  (two  equation  models,  eddy  viscosity 
models)  will  have  a  very  difficult  time  modeling  rotating  turbulence.  Mathematically,  it  can  also  be  shown 
that  the  turbulent  kinetic  energy  alone  is  not  a  sufficient  descriptor  of  the  turbulence  when  rotation  is 
present.  The  more  complex  Reynolds  stress  transport  (RST)  models  appear  to  be  the  simplest  turbulence 
models  that  are  physically  rich  enough  to  be  able  to  capture  rotation.  These  models  contain  explicit  and 
exact  Coriolis  terms,  as  well  as  some  rotational  terms  which  need  to  be  modeled. 

The  turbulent  potential  model  is  essentially  a  reduced  (less  expensive)  RST  model.  It  is  expected  therefore 
that  this  model  will  also  be  able  to  capture  rotating  turbulent  flows  reasonably  well.  However,  the  slight 
disadvantage  of  the  potential  model  approach  is  that  the  Production  and  Coriolis  terms  are  only  exact  in 
flows  with  a  single  Cartesian  inhomogeneous  direction.  This  is  apparently  the  price  to  be  paid  for  the 
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reduced  cost.  In  flows  with  complex  inhomogeneity  in  the  turbulence,  which  frequently  is  the  case  in 
rotating  flows,  some  additional  modeling  is  required  over  pure  RST  models. 

The  modeling  of  the  production  and  Coriolis  terms  is  described  in  detail  below.  The  key  principal  in  this 
modeling  effort  is  the  idea  that  the  equations  should  produce  the  same  answer  if  you  are  in  a  fixed  reference 
frame  that  sees  a  rotating  mean  flow  or  if  you  are  in  a  rotating  reference  frame  that  moves  with  the  mean 
flow  and  therefore  has  zero  mean  flow. 


2.2  Production 

The  production  terms  and  have  to  be  modeled  in  rotating  flows.  In  order  to  model  the  production 
terms,  we  look  at  the  form  of  the  terms  in  flows  with  only  one  direction  of  inhomogeneity,  and  then  try  to 
generalize  the  terms  for  all  flows.  We  first  split  the  tensor  Py  into  two  parts 


P.  =  P?  +  P 

y  y  y 


w 


(2.1) 


Here 


Pjj  “  R&Sjk  RjkSik 

P*w  =  -RikWjk  -  Rjkwik 

For  flows  lying  in  the  xrx2  plane  with  inhomogeneity  in  only  the  x2  direction,  we  therefore  write 


(2.2) 


P,S2=-fl+-^ 


R 


p  C  pw  — 

JV22°12  r12  “ 


f  P  ^ 

\  _  J^ji 


22  J 


R 


r22w12 


22  J 


^22  “  ^^12^12 


PW=2R  W 

r22  ^'*'‘12  vv12 


(2.3) 


Using  the  fact  that  for  these  kind  of  flows  (1 .9)  holds,  we  can  then  say 


=  (1 .  -  q(V>  <t>>  k))4>s3  PV3  =  q(¥?  k)(()a): 


P*  =  V3S3  P*W=-V3®3 


(2.4) 


1 

Where  q  is  a  scalar  meant  to  represent  the  ratio  — 


j  Rn 

R 


A 


,  which  vanishes  at  isotropy,  s  is  a  vector 


22  J 


representing  strain  and  0)  is  a  vector  representing  intrinsic  vorticity.  Now  it  is  easy  to  see  from  (2.3)  and 
(2.4)  that  the  vectors  have  to  obey  the  following  relationship  with  W  and  S 


6)3  =  -2W,2 

s3  =  — 2SI2 


(2.5) 


For  3D  flows  with  inhomogeneity  in  only  the  x2  direction,  an  additional  set  of  relationships  arise 


(bj  =  2W32 

Sj  =  2S32 


(2.6) 
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We  now  try  to  define  0)  and  s  more  generally  so  that  (2.5)  and  (2.6)  can  be  satisfied  for  3D  flows  with 
inhomogeneity  in  only  the  x2  direction.  The  definition  of  (5)  is  quite  straightforward 

<$)=Vxv  +  2«  (2.7) 

Which  is  the  same  as  the  intrinsic  vorticity  defined  by  Speziale  (1989).  However,  the  relationship  between 
s  and  the  gradients  of  velocity  is  not  algebraic,  and  the  following  equation  is  solved 

V2  (<{>s)  =  -2Vx  V  •  (4>S)  (2.8) 


This  equation  ensures  that  s  satisfies  the  relations  (2.5)  and  (2.6)  for  the  specific  case  of  flows  with  only 
one  direction  of  inhomogeneity.  However,  for  general  flows,  we  need  to  solve  for  the  inverse  Laplacian  of 
the  source  term  defined  on  the  right  hand  side  of  (2.8).  For  wall  bounded  flows,  (])  goes  to  zero  at  the  wall, 
and  therefore  the  following  boundary  condition  would  be  appropriate 

<Hwa„=0- 

Finally  we  can  define  the  production  terms  for  the  turbulence  potentials 

p*  =  VA  P«w  =  -Vk®k  P*  =  VkSk  -  Vk®k 

Fj  =  (1  -  q(V>  <|>,  k))<j)s,  Pvw  =  q(y, <>,  k)(J)6)i 
PVi  =  (1  -  q(y>  k))<f>si  +  q(\|f,  <j>,  k)<^ 

The  scalar  function  q  has  been  defined  in  the  section  defining  the  model. 

2.3  Decaying  Homogenous  Isotropic  Turbulence  in  a  Rotating  Frame 

The  prediction  of  decaying  isotropic  turbulence  in  a  rotating  frame  can  be  used  to  calibrate  the  k-e 
equation,  more  specifically  Ce2.  Experimentally,  rotating  homogenous  isotropic  turbulence  was  produced 
by  Wigeland  and  Nagib  (1978)  by  setting  a  uniform  flow  into  solid  body  rotation  as  it  passed  through  a 
honeycomb  and  a  rotating  turbulence-  generating  grid.  It  has  been  known  that  the  dissipation  rate  decreases 
rapidly  due  to  a  disruption  of  the  energy  transfer  from  larger  to  the  smaller  scales  Speziale  (1989).  C£2  also 
has  a  dependence  on  the  Reynolds  number  as  seen  by  Chasnov  (1997).  The  effect  of  the  Reynolds  number 
and  rotation  can  be  seen  in  (1.14),  so  that  an  increase  in  rotation  causes  an  increase  in  Cr  while  an  increase 
in  Reynolds  number  makes  Ce2  tend  to  go  closer  to  1.83.  Ce2  has  a  lower  limit  of  1.664  for  an  initial 
condition  with  a  k2  energy  spectrum  at  small  wave  numbers. 

Figure  1  shows  the  results  of  the  prediction  of  isotropic  decay  for  3  cases  tried  out  by  Wigeland  and  Nagib 
(1978).  The  x-axis  is  nondimensional  time  Ut/M,  where  U  is  the  mean  streamwise  velocity  of  the  air,  M  is 
the  size  of  the  mesh  which  generates  the  turbulence,  and  t  is  the  time.  The  y-axis  is  the  nondimensional 
reciprocal  of  turbulent  kinetic  energy  kO/k,  where  kO  is  the  initial  k.  The  reciprocal  is  taken  to  simply 
resolve  the  k  more  effectively  for  different  cases.  Ret  =  ko2A)£o>  the  initial  turbulent  Reynolds  number, 
while  Rot=8o/ftko,  the  initial  turbulent  Rossby  number. 

As  can  be  seen  from  the  results,  the  predictions  match  quite  accurately  with  the  model,  and  therefore  we 
can  say  that  Ce2  has  been  calibrated  correctly  within  the  range  of  Reynolds  and  Rossby  numbers  displayed 
in  the  figures.  The  Reynolds  number  used  in  later  flows  are  about  10  times  larger  than  what  has  been  used 


(2.9) 


(2.10) 
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in  these  cases.  However,  it  was  seen  that  isotropic  decay  is  more  sensitive  to  lower  Reynolds  numbers,  and 
is  independent  of  Re,  for  high  values  of  Re,. 

The  effect  of  Rossby  number  on  the  evolution  of  isotropic  turbulence  at  high  Reynolds  number  can  be  very 
clearly  seen  in  the  last  plot  in  Figure  1.  The  model  is  compared  to  a  high  Reynolds  number  DNS  presented 
in  Speziale  (1989),  and  depicts  the  effect  of  the  rotation  dependent  factor  on  the  decay  of  turbulence.  The 
initial  turbulent  Reynolds  number  is  390  and  the  initial  turbulent  Rossby  number  is  0.01 1.  The  model  has 
been  tried  with  and  without  the  rotation  dependent  factor  Cr  (equation  1.14).  Clearly,  the  presence  of  Cr 
decreases  the  dissipation  rate  to  the  right  amount  so  that  the  model  agrees  well  with  the  experiment. 


Case  A  Cass B 


Case  C 


Figure  1.  Cases  A  B  and  C  of  Wigeland  and  Nagib’s  (1978)  experiment.  The  symbols  denote  the 
experiment  while  lines  denote  model  predictions.  Final  figure  is  the  DNS  results  of  Speziale  at  Ro=0.01 1. 


2.4  Rotating  Homogenous  Planar  Shear  Flow 


In  this  test  case,  initially  homogenous  isotropic  turbulence  sheared  in  a  plane  at  a  constant  rate  and  a  fixed 
amount  of  frame  rotation  is  simultaneously  superimposed.  The  rotation  is  around  an  axis  perpendicular  to 
the  plane  of  the  shear.  The  schematic  is  shown  in  Figure  2. 
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Figure  2.  Homogenous  shear  flow  (black  arrows)  superimposed  with  frame  rotation  (block  arrow), 
acting  on  initially  isotropic  turbulence.  In  this  figure  the  shear  and  the  rotation  have  the  same  sign. 


All  the  equations  come  into  play,  and  only  terms  involving  the  gradient  of  turbulence  potentials  are  equal  to 
zero  because  of  the  homogeniety.  The  shear  is  denoted  by  S  =  & .  The  rotation  rate  by  Q.  Initially, 

<|>=0.66k,  \|/3  =0,  and  £o/k0S  =  0.25.  Zofk 0S  &  Q/S  are  the  only  initial  condition  on  which  the  solution 
depends,  as  shown  by  Speziale  &  Mac  Giolla  Mhuiris  (1989). 

The  prediction  of  homogenous  shear  flow  without  frame  rotation  is  compared  to  DNS  data  by  Matsumoto, 
Nagano  and  Tsuji  (1991).  From  Figure  3  we  can  see  that  the  model  predicts  a  lower  turbulent  kinetic 


Figure  3.  Evolution  of  turbulent  kinetic  energy  with  time  for  homogenous  shear  flow  with  no  rotation. 
Symbols  denote  DNS  data.  Line  denotes  model  prediction 
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energy  which  seems  to  be  going  towards  a  stable  value  for  large  times.  The  initial  evolution  does  match 
however,  and  therefore  the  value  of  0.6  (at  isotropy)  chosen  for  Csp2  is  correct. 

From  Figure  4  we  can  see  that  the  model  agrees  reasonably  well  with  the  LES  for  homogenous  shear  flow 
with  frame  rotation,  and  predicts  stable  solution  for  Q/S— -0.5.  The  turbulent  kinetic  energy  is  unstable  for  a 
range  of  Q/S  between  -0.1  and  0.65,  which  is  higher  than  the  range  of  0.<  Q/S<0.5  predicted  by  linear 
stability  theory  and  RDT  (Lezius  and  Johnston,  1976  and  Bertoglio,  1982). 


Figure  4.  Effect  of  homogenous  shear  and  rotation  on  initially  isotropic  turbulence.  Symbols  denote  LES 
data  of  Bardina  et  al  (1983)  and  lines  denote  model  predictions.  The  symbol  O  is  for  Q/S=0.5,  A  is  for 
Q/S=-0.5 


Thus  the  model  does  not  predict  the  homogenous  shear  flow  very  well  because  it  predicts  a  slow  evolution 
rate  for  no  frame  rotation  rate.  This  is  of  some  concern  because  this  evolution  rate  is  shown  to  be  predicted 
correctly  for  the  relatively  simple  k-E  model  (Speziale,  1993). 


2.5  Rotating  Spanwise  Channel  Flow 


The  schematic  for  this  case  is  shown  in  Figure  5.  A  fully  developed  flow  between  two  infinitely  wide  and 
deep  plates  is  rotated  around  an  axis  perpendicular  to  the  plane  of  the  flow.  The  predictions  for  this  case  are 
compared  to  DNS  by  Kristoffersen  and  Andersson  (1993).  The  DNS  has  a  turbulent  Reynolds  number 


Ret  =  UTh/t)  of  194.  Here  h  is  the  channel  half  width  and  uT  = 


dv 

*  U — 

V  dy 


wall 


is  the  turbulent  velocity  scale. 


The  pressure  gradient  is  equal  to  1.  The  Rossby  number  Ro  for  this 


case  is  defined  as  2h|fi|/Ub ,  where 


Ub  is  the  bulk  mean  velocity  of  the  flow  in  the  channel. 


The  primary  issue  in  modeling  the  equations  for  turbulent  potentials  for  the  case  of  channel  flow  was  to  get 
the  correct  wall  asymptotics.  In  the  equation  for  <j),  if  Cp2  and  C^2  tend  to  some  arbitrary  constant  other 

than  2.0  near  the  wall,  then  the  terms  involving  the  strain  vector  and  the  intrinsic  vorticity  vector  vary  like 
0(y3)  near  the  wall  while  the  other  terms  vary  like  0(y5)  close  to  the  wall.  This  leads  to  difficulty  in 

integrating  the  equation  up  to  the  wall.  In  order  to  avoid  this  we  make  Cp2  and  Cp2  tend  to  2.0  near  the 


9 


Suction  side 


Figure  5.  Fully  developed  spanwise  rotating  channel  flow.  The  Block  arrow  denotes  the  direction  of 
positive  rotation.  The  axis  of  rotation  is  at  a  very  large  distance  compared  to  the  channel  width. 

wall  in  such  a  way  such  that  an  extra  factor  of  (1  -  a)  now  multiplies  the  terms  involving  the  intrinsic 
vorticity  vector  and  the  strain  vector.  Since  (1  —  OC)  varies  near  the  wall  as  0(y2),  therefore  the  wall 
asymptotics  of  the  <j)  equation  is  correct,  and  we  can  integrate  the  equation  up  to  the  wall. 

From  Figure  6  we  can  see  that  for  zero  rotation,  the  mean  velocity  field  is  well  predicted.  Figure  7  shows 
that  the  shape  of  the  velocity  profile  becomes  assymmetric  for  Ro  =  0.15,  though  the  slope  of  the  velocity 
profile  predicted  on  the  pressure  side  is  not  correct.  Figure  8  shows  that  the  slope  is  predicted  correctly  for 
a  higher  rotation  rate  of  Ro  =  0.5.  Therefore  we  can  say  that  the  turbulent  potential  model  is  sensitive  to 
spanwise  rotation  of  a  fully  developed  channel,  and  gives  a  good  prediction  of  the  velocity  profile 


y/h 


Figure  6.  Mean  velocity  profile  for  Ret  =  194.  No  rotation  is  present.  Mean  velocity  has  been 
nondimensionalized  with  respect  to  UT .  Symbols  denote  DNS  data  while  line  denotes  model  prediction. 
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2.6  Isolated  turbulent  trailing  vortex 

A  turbulent  trailing  vortex  is  created  by  the  wingtip  of  a  plane.  Since  the  trailing  vortex  is  a  case  for  which 
few  experimental  or  numerical  results  are  present,  therefore  we  make  some  assumptions  about  the  vortex 
which  are  the  same  as  those  made  by  Wallin  and  Girimaji  (2000).  Cylindrical  coordinates  are  used  to 
represent  the  components  of  the  velocity  field.  Initially  both  mean  axial  and  azimuthal  velocities  are 
present  in  the  vortex.  After  some  time,  assuming  that  the  axial  velocity  diffuses  due  to  the  turbulence,  only 
the  azimuthal  velocity  will  be  present.  Consequentially  the  flow  is  homogenous  in  the  z  direction.  The 

mean  azimuthal  velocity  (i.e.  V0)  gets  affected  only  by  the  component  \|/, .  In  order  that  the  vortex  grows 
at  the  right  rate,  \| /,  should  have  the  correct  magnitude  and  radial  profile. 

A  major  difference  between  this  case  and  the  previous  cases  is  that  the  turbulence  is  inhomogeneous  in  2 
directions,  and  therefore  there  is  no  one-to-one  correspondence  between  the  turbulent  potentials  and 
selected  components  of  Reynolds  stresses.  Therefore,  this  is  the  first  case  where  we  will  numerically  solve 
(2.8)  and  find  out  the  component  of  the  strain  vector  sh  Since  this  is  not  a  wall  bound  flow  therefore  we 
cannot  take  (2.9)  as  a  boundary  condition.  However,  we  can  assume  that  at  very  large  radius  the  turbulent 
quantities  are  inhomogeneous  only  in  the  r  direction,  i.e.  the  cylindrical  and  Cartesian  coordinates  start 
becoming  similar  at  large  radius.  We  can  also  assume  close  to  isotropic  conditions  at  large  radius. 
Therefore  we  can  take  our  boundary  conditions  as 


Along  with  (2.1 1)  we  take  zero-derivative  boundary  conditions  for  the  turbulent  potentials.  The  results  for 
this  problem  are  presented  in  terms  of  the  circulation  T  =  27H’V0  .  Potential  flow  is  present  at  points  far 
away  from  the  vortex  core  because  of  the  absence  of  any  significant  stresses.  Therefore,  the  circulation  will 
be  a  constant  T0  at  the  boundary.  Hence,  we  take  the  zero-derivative  boundary  condition  for  the 

circulation  as  well.  LIDAR  (Coherent  Doppler  Laser  Radar)  measurements  (Campbell  et  al.,  1996)  of  an 
aircraft  trailing  vortex  at  Memphis  airport  have  been  taken  for  initial  conditions  and  for  comparing  the  final 
conditions  with  the  results  predicted  by  the  model.  The  initial  conditions  that  have  been  assumed  for  the 
turbulent  quantities  are  taken  according  to  the  suggested  values  given  by  Wallin  and  Girimaji  (2000).  All 
the  initial  conditions  are  nondimensionalized  with  respect  to  F0=350  m  /s  and  a  critical  radius  Rc=L8  m. 

They  are  ko=10'5  and  e0=10'8.  We  also  assume  (t)o=0.66ko  and\|/j=0  because  of  initially  isotropic 

conditions.  Turbulent  Reynolds  number  Re^OX  106.  The  LIDAR  measurement  for  the  final  circulation 
profile  has  been  taken  at  a  time  t  =  55  s  (which  is  equivalent  to  about  6000  non-dimensional  time  steps). 
The  final  result  for  this  case  can  be  seen  in  Figure  9.  The  circulation  is  falling  at  around  log(r)~0.3,  as 
apparent  in  the  experiment.  However,  the  rate  of  decrease  is  much  slower  than  the  measured  values.  Also, 
there  seems  to  be  an  overshoot  of  the  circulation,  which  is  physically  possible  (Govindaraju  and  Saffman, 
1971),  it  is  not  true  for  the  experimental  results.  The  circulation  at  around  log(r)=0.1  seems  to  have 
increased  instead  of  decreasing. 

Thus  the  results  for  the  turbulent  trailing  vortex  are  not  highly  satisfactory,  as  the  growth  rate  of  the  vortex 
seems  to  be  very  low.  However,  the  model  prediction  does  have  the  feature  of  having  an  unchanged  core 
region.  This  is  significantly  better  than  the  k-£  model  which  gives  a  rapid  expansion  of  the  core  region 
(Wallin  and  Girimaji,  2000). 
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1.2 


Figure  9.  Circulation  profiles  of  a  trailing  vortex.  Symbol  *  indicates  initial  circulation  profile  for  both 
model  and  LIDAR  measurement.  Solid  line  gives  final  circulation  profile  predicted  by  model  at  t=55  s. 
Symbol  o  indicates  final  circulation  profile  measured  by  LIDAR 

2.7  Swirling  Pipe  Flow 

The  swirling  pipe  flow  is  the  most  difficult  case  to  model  because  the  one-to-one  relationship  between  the 
Reynolds  stress  and  the  turbulence  potentials  is  lost.  Figure  10  shows  the  schematic  diagram  for  this  case. 
The  z,  r,  0  directions  correspond  to  1,  2  and  3  directions  respectively. 


Figure  10.  Schematic  of  swirling  pipe  flow.  The  velocity  field  is  shown  relative  to  a  stationary  frame  of 
reference. 
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The  model  has  been  compared  to  the  experiment  by  Imao  and  Itoh  (1996).  The  calculations  have  been  done 
with  respect  to  a  non-rotating  frame,  so  that  Coriolis  terms  are  not  present  in  the  calculations.  The  radius  of 
the  pipe  is  R.  The  turbulent  Reynolds  number  u,RA)  for  this  case  is  572.  The  rotation  rate  is  defined  as 
vj  /v  I  ,  and  the  three  rotation  rates  tried  out  are  0,  0.5  and  1.0.  The  non-dimensional  pressure 

“I  wall  2 1  mean 

gradients  applied  for  the  three  cases  are  1.0,  0.8  and  0.6  respectively. 

The  extra  term  which  contributes  to  the  model  for  this  case  is  the  vector  A,  defined  in  (1.14).  This  term 
does  not  contribute  to  the  channel  flow  or  2-D  flows.  This  is  because  the  term  consists  of  vectors  like 

S  X  \J /  and  (6)-  s)  X\j/  which  is  zero  if  the  vectors  G)  and  s  are  parallel  to  .  However,  for  a  3D  flow 
this  term  can  affect  all  the  components  of  \J f . 

Figure  11  compares  the  axial  and  azimuthal  velocity  profiles  for  N  equal  to  0.,  0.5  and  1.0  respectively. 
Although  the  agreement  of  the  model  predictions  with  the  experimental  data  is  not  very  good,  it  does 
follow  the  trends  shown  by  the  data.  As  we  go  to  higher  rotation  rates,  the  centerline  axial  velocity 
increases.  Also,  as  rotation  rate  increases,  the  deficit  between  the  azimuthal  velocity  profile  predicted  by 
the  model  and  the  laminar  (and  linear)  azimuthal  velocity  profile  increases.  The  overshoot  of  the  azimuthal 
velocity  profile  near  the  wall  decreases  as  the  rotation  rate  is  increased. 
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Figure  1 1.  (a)  Axial  velocity  profile  for  N=0,  (b)  N=0.5,  (c)  N=1.0.  (b)  Axial  and  azimuthal  velocity  profiles 
for  N  =  0.5.  (c)  Axial  and  azimuthal  velocity  profiles  for  N  =  1.0.  Symbols  +  indicate  experimental  results 
for  axial  velocity.  Symbols  o  indicate  experimental  results  for  azimuthal  velocity.  Bold  lines  indicate  model 
predictions.  Dashed  lines  depict  azimuthal  velocity  profile  for  laminar  flow. 
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Nondimensional  spanwise  velocity  Nondimensional  axial  velocity 


Figure  12.  Schematic  of  the  streamwise  rotating  channel  flow 


Figure  13.  (a)  Axial  mean  velocity  at  Ro  =  3.  (b)  Axial  mean  velocity  at  Ro  =  10.  (c)  Spanwise  mean  velocity 
at  Ro  =  3.  (d)  Spanwise  mean  velocity  at  Ro  =10.  Symbols  denote  DNS.  Line  denotes  model  prediction. 


The  streamwise  rotating  channel  flow  consists  of  fully  developed  flow  between  two  channels  in  the  Xi 
direction  which  is  then  rotated  around  an  axis  along  the  xl  direction.  As  shown  in  Figure  12  both  axial  and 
spanwise  velocities  get  generated  in  this  flow.  There  is  only  one  direction  of  inhomogeneity  (i.e.  x2).  Thus 
we  have  a  one-to-one  correspondence  between  some  of  the  Reynolds  stresses  and  the  turbulent  potentials. 
This  case  has  been  tried  out  along  with  the  swirling  pipe  flow  in  order  to  correctly  calibrate  the  vector  A  in 
the  equation  for  \| /  .  The  vector  A  will  appear  in  this  case  because  the  flow  is  3D.  The  vector  A  is  not  seen 
in  the  equation  for  turbulent  shear  stresses  in  RANS  models,  although  a  one-to-one  correspondence  exists 
between  the  vector  \|f  and  the  shear  stresses  R12,  R23.  This  may  be  because  the  exact  equations  for  the 
turbulent  shear  stresses  involve  terms  containing  factors  like  (Rn  -  R22)  and  (R33  -  R22),  which  cannot  be 
represented  by  only  and  k.  The  results  of  the  model  predictions  have  been  compared  to  DNS  by  Oberlack, 
Cabot  and  Rogers  (1998).  The  Reynolds  number  for  the  flow  is  given  by  Ret  =  hux  / 1)  =  1 80  where  all 
the  quantities  have  the  same  definition  as  for  the  case  of  spanwise  rotating  channel  flow.  The  Rossby 
number  is  given  by  Ro  =  2Qh/ux.  The  DNS  data  is  there  for  Ro=3.0  and  Ro=  10.0,  therefore  these  are 

the  rotation  rates  against  which  the  model  is  compared.  From  (a)  Axial  mean  velocity  at  Ro  =  3.  and 
Error!  Reference  source  not  found,  we  can  see  that  the  streamwise  velocity  is  predicted  quite  well  by  the 
model.  For  higher  rotation  rates  the  model  seems  to  be  having  a  tendency  to  give  a  more  parabolic  profile 
at  the  center,  whereas  the  DNS  seems  to  be  giving  a  flatter  profile.  From  Figure  13  we  see  that  the 
spanwise  velocity  is  predicted  well  for  low  rotation  rates,  but  deviates  appreciably  from  the  DNS  at  high 
rotation  rates. 

2.9  Conclusions 

•  The  model  correctly  predicts  that  the  dissipation  rate  decreases  with  an  increase  in  rotation  rate. 

•  The  model  is  frame  consistent  with  frame  invariant  flow-dependent  vectors  mathematically  defined. 
The  frame  rotation  rate  appears  in  the  equations  only  through  the  intrinsic  vorticity  vector. 

•  The  mode!  gives  the  wrong  evolution  rate  for  the  turbulent  kinetic  energy  for  initially  isotropic 
turbulence  subjected  to  plane  shear. 

•  2D  rotating  channel  flow  has  been  well  predicted. 

•  Model  predictions  do  not  give  the  required  growth  rate  for  the  trailing  vortex.  However,  the  results  are 
better  than  the  k-e  model. 

•  A  vector  A  is  present  in  the  model  which  is  constructed  in  such  a  way  that  it  contributes  only  for  3D 
flows.  Both  swirling  pipe  flow  and  streamwise  rotating  channel  flow  have  been  predicted  reasonably 
well. 

The  turbulent  potential  model  has  been  constructed  in  such  a  way  that  it  is  sensitive  to  frame  rotation  as 
well  as  streamline  curvature.  At  the  same  time,  it  is  also  frame  consistent.  In  order  to  construct  a  frame 
consistent  model,  we  first  identified  the  form  of  the  transport  equation  for  the  turbulent  potentials,  and  then 
constructed  frame  invariant  flow  dependent  vectors  (i.e.  the  strain  and  intrinsic  vorticity  vectors).  Thus  we 
ensure  that  any  dependence  of  the  model  on  the  frame  rotation  rate  has  to  be  through  the  intrinsic  vorticity 
vector.  The  similarity  of  the  model  with  RANS  models  for  flows  with  only  one  direction  of  inhomogeniety 
is  then  used  in  order  to  find  out  the  form  of  the  various  budget  terms  in  the  model  (i.e.  pressure  strain, 
dissipation,  turbulent  transport  etc.).  This  model  is  then  extended  for  3D  flows,  and  an  extra  term  has  to  be 
added  which  contributes  to  only  3D  flows. 

The  results  show  that  the  model  does  a  good  job  of  predicting  the  flows  in  the  axial  direction  for 
swirling/streamwise  rotating  flows.  The  model  also  predicts  correctly  the  streamwise  velocity  profile  for 
the  spanwise  rotating  channel  flow.  The  spanwise/azimuthal  velocity  is  not  well  predicted  for  either 
swirling  flows  or  the  trailing  vortex  flow,  however  the  spanwise/azimuthal  velocities  are  much  smaller  than 
the  axial  velocity  and  therefore  it  is  not  as  important  as  the  axial  velocity.  The  wrong  evolution  rate  for  the 
non-rotating  homogenous  shear  flow  is  of  some  concern. 
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Predicting  rotating  turbulent  flows  is  a  very  difficult  task.  The  corrections  to  the  model  to  account  for 
rotation  are  not  as  elegant  as  initially  desired.  Nevertheless,  a  great  deal  was  learned  in  this  project 
concerning  rotating  turbulent  flows,  their  physics,  and  how  to  correctly  model  rotation. 


3.  Transition 


3.1  Introduction 

Transition  to  turbulence  is  a  stochastic  process  which  is  only  partially  understood.  Depending  to  a  large 
extent  on  the  characteristics  of  the  initial  disturbances  a  laminar  boundary  layer  flow  will  experience 
different  transition  processes  until  it  finally  develops  into  fully  developed  turbulence.  Typically,  two  types 
of  transition  scenarios  are  considered  -  natural  transition  and  bypass  transition. 

If  the  magnitude  of  the  disturbance  is  small  enough,  the  excited  instabilities  will  grow  like  typical  linear 
theories  predict.  Once  the  instabilities  are  large  enough,  spanwise  fluctuations  will  be  triggered  and  then 
develop  to  the  same  scale  as  the  streamwise  fluctuation  waves.  Subsequently,  non-linear  effects  take  over 
and  transition  will  take  a  route  of  exponential  growth  of  instability  waves,  which  is  known  as  “secondary 
instability”  and  is  a  result  of  the  interaction  between  the  inertial  scale  disturbances  and  local 
heterogeneities.  Following  this  stage,  the  process  is  rather  a  breakdown  than  a  growth.  Instabilities 
ultimately  burst  forth  into  patches  of  irregular  motion  which  are  termed  turbulence  spots  (Emmons,  1951). 
These  intermittent  spots  are  interspersed  by  laminar  flow  but  they  develop  rapidly  and  eventually  coalesce 
into  fully  turbulent  flow.  A  transition  process  evolving  through  the  stages  described  above  is  termed 
“natural  transition”. 

Alternatively,  in  some  circumstances,  the  disturbances  can  be  at  the  inertial  scale  initially  and  the  non¬ 
linear  effects  immediately  come  in  to  play.  The  linear  growth  process  is  bypassed  and  secondary  instability 
and  turbulence  spot  occurs  very  quickly  and  the  flow  becomes  turbulent  in  a  short  distance.  Obviously, 
such  a  transition  mechanism  is  quite  different  from  natural  transition.  It  is  termed  “bypass  transition”  and  is 
stochastic  in  nature.  In  addition  to  the  difference  in  stages  experienced  by  the  flow,  the  characteristics  of 
the  flow  patterns  in  their  common  stages  can  be  different.  For  example,  for  a  transitional  boundary  flow, 
the  turbulent  spots  in  the  bypass  transition  are  generated  by  free-stream  eddies  while  the  classical  Emmons 
type  spots  in  the  natural  transition  are  usually  stimulated  by  instabilities  near  the  wall  (Jacobs  &  Durbin, 
2001). 

Due  to  the  existence  of  apparently  different  transition  mechanisms,  the  attempts  to  develop  a  universal 
methodology  to  model  transition  to  turbulence  are  often  met  with  skepticism.  Nevertheless,  for  most 
engineering  applications,  we  only  require  the  model  to  produce  an  accurate  prediction  of  the  mean  flow 
behavior  in  the  transition  region.  This  allows  us  to  neglect  the  mechanism  of  instability  evolution  in  the 
transition  and  only  take  care  of  quantities  that  contribute  significantly  to  the  mean  flow  variables.  There  is 
no  doubt  that  we  can  always  find  some  general  traits  of  flows  in  the  transitional  zone,  which  distinguish 
them  from  either  laminar  or  full  turbulent  flows.  Basically,  the  transitional  zone  is  a  region  for  the 
fluctuations  of  flow  variables  to  develop  from  small  scales  into  the  full  span  of  turbulent  scales.  These 
intensities  of  the  fluctuating  quantities  are  large  enough  to  measure  the  “random”  behavior  of  the  flow  and 
they  are  small  enough  to  keep  the  stochastic  part  from  affecting  the  mean  flow  significantly  until  the  very 
end  of  the  transition.  The  methodology  of  this  work  recognizes  the  fact  that  the  contribution  to  the  mean 
flow  of  the  small  amplitude  fluctuations  in  transitioning  flow  can  be  taken  into  account  by  certain  terms  in 
the  turbulence  model. 

3.2  Review  of  Transition  Prediction  Approaches 

Predicting  the  onset  of  turbulent  flow  is  a  critical  component  of  many  engineering  and  environmental 
flows.  The  characteristics  of  laminar  and  turbulent  boundary  layers  are  so  different  that  the  precise  location 
of  this  relatively  abrupt  transition  can  have  a  profound  influence  on  the  overall  drag,  heat  transfer,  and 
performance  properties  of  devices  that  operate  in  the  transitional  regime.  The  prediction  of  boundary  layer 
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transition  is  complicated  by  the  fact  that  it  does  not  correspond  very  directly  to  the  onset  of  instability. 
Stability  analysis  for  boundary  layers  is  well  developed  and  very  predictive  of  the  behavior  of  small 
disturbances,  however  the  instabilities  go  through  a  series  of  complex  non-linear  and  three-dimensional 
processes  before  turbulence  itself  actually  develops.  In  addition,  stability  analysis  is  less  helpful  with  the 
prediction  of  by-pass  transition  where  external  free-stream  disturbances  by-pass  the  classic  instability 
mechanisms  and  initiate  the  non-linear  three-dimensional  transition  process  directly. 

Traditional  methods  for  predicting  transition  rely  on  correlations.  For  example,  van  Driest  &  Blumer 
(1963)  suggest  the  implicit  correlation: 

1690  Re'1,,2  =  0.3 12(m  +  0. 1 1)"'528  +  4.8J929  Re]/2  T2  (3.1) 


where  Re^  tr  is  the  transition  Reynolds  number  based  on  the  local  free-stream  velocity  U ,  the  parameter  m 
is  related  to  the  dimensionless  pressure  gradient  (m=0  for  zero  pressure-gradient),  d99  is  the  local  99% 

boundary  layer  thickness,  and  T  =  yf^k  !U  is  the  local  free-stream  turbulence  intensity.  This  is  only 

one  of  many  proposed  correlations.  While  they  appear  relatively  simple,  they  are  actually  quite  awkward 
to  implement  in  a  general  purpose  CFD  code.  In  a  general  code,  correlations  require  each  point  on  the 
boundary  to  determine  non-local  quantities,  such  as  free-stream  values,  distance  downstream  from  the 
leading  edge  and  boundary  layer  thickness.  Unambiguous  definitions  for  such  quantities  in  a  general 
situation  are  very  difficult  to  formulate.  Many  CFD  codes  simply  resort  to  requiring  the  user  to  specify  the 
transition  location.  Once  the  transition  location  has  been  determined,  there  is  the  further  difficulty  of 
determining  how  the  turbulence  model  should  be  prompted  to  become  active.  Solutions  to  this  problem 
range  from  specifying  large  turbulent  kinetic  energy  at  the  wall  where  transition  is  expected  to  occur  to 
specifying  local  large  source  terms  in  the  turbulence  evolution  equations  (artificial  production).  None  of 
these  methods  of  ‘tripping’  the  turbulence  model  are  particularly  reflective  of  the  actual  physical  transition 
process,  and  add  further  ambiguity  to  the  already  uncertain  transition  location.  In  addition,  ‘tripping’  of 
this  sort  is  largely  code  dependent,  so  that  the  same  theoretical  model  can  produce  different  results  in 
different  codes,  or  even  within  the  same  code  using  different  meshes. 

An  alternative  approach  is  to  use  the  turbulence  model  itself  to  predict  the  transition  location.  This  is  a  very 
natural  approach  in  by-pass  transition,  since  the  model  is  effectively  on  in  the  free-stream  anyway  to 
predict  the  free-stream  turbulence.  A  number  of  studies  on  boundary  layer  by-pass  transition  prediction 
using  low  Reynolds  number  k/E  models  have  been  performed,  and  Savill  (1990  &  1996)  are  good 
reviews  of  how  various  flavors  of  the  k  /  €  models  perform.  The  overall  conclusion  is  that  none  performs 
well  for  all  the  flows  considered.  More  recent  work  using  a  k  I  £  intermittency  model  (Suzen  &  Huang, 
2001)  has  shown  more  success.  However,  none  of  these  models  attempts  to  predict  natural  transition  or 
relaminarization. 


3.3  Motivation 

While  the  idea  of  using  RANS  models  to  predict  transition  may  have  originally  been  motivated  by  practical 
considerations  such  as  ease  of  implementation,  it  also  has  a  solid  theoretical  justification  that  has  not  been 
discussed  previously.  The  derivation  of  the  Reynolds  Averaged  Navier-Stokes  (RANS)  and  accompanying 
Reynolds  Stress  Transport  (RST)  equations  are  simply  a  mathematical  reformulation  of  the  governing 
Navier-Stokes  equations.  There  are  no  physical  assumptions  or  restrictions  on  the  range  of  applicability  of 
these  equations,  as  presented  below. 

•^  +  V  •  (uu)  =  -Vp  +  V  •  vVu  -  V  •  R  (3.2a) 
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turbulent  transport,  and  dissipation  are  only  expected  to  influence  the  very  last  stages  of  transition  (the 
overshoot  recovery)  and  should  have  little  affect  on  the  critical  transition  location.  Most  models  for  the 
rapid  pressure  strain  usually  only  capture  rapid  shear  correctly  because  this  is  the  important  case  for 
boundary  layer  and  free-shear  flows.  Reynolds  and  Kassinos  (1995)  have  a  modeling  framework  that  can 
capture  any  rapidly  distorting  flow. 

Due  to  the  complexity  of  solving  the  Reynolds  stress  transport  equations,  two-equation  models  such  as 
k  /  £  are  very  popular.  The  k-equation  is  derived  by  taking  one  half  of  the  trace  of  Eqn  3.2. 


dK_ 

dt 


1 


+  u -VK  =  — u'u 
2  J 


dx. 


rd»-y 


\?xj  j 


+  vV2K 


(3-4) 


Two  very  important  changes  take  place  when  only  the  k-equation  is  solved.  First,  pressure  strain 
disappears  entirely.  This  is  because  this  term  only  redistributes  energy  between  various  Reynolds  stress 
components,  it  does  not  influence  the  overall  energy  content.  Unlike  spatial  redistribution  (transport  terms) 
this  intercomponent  redistribution  is  important,  particularly  in  the  boundary  layer.  It  takes  energy  out  of 
the  amplified  streamwise  fluctuations  and  transfers  it  to  the  normal  and  spanwise  fluctuations,  damping  the 
streamwise  fluctuation  growth  and  feeding  the  normal  and  spanwise  velocity  fluctuations.  Two  equation 
models  (like  k/£)  are  fundamentally  incapable  of  capturing  physical  processes  that  involve 
intercomponent  energy  redistribution.  In  addition,  the  production  term  is  no  longer  exact.  The  Reynolds 
stress  tensor  in  this  term  must  now  be  modeled.  The  classic  Boussinesq  eddy  viscosity  model  is 

R  =  |/il-Kr(Vu  +  Vur)  (3.5) 


C  K 2 

Where  the  eddy  viscosity  is  given  by  yf  =  - More  complicated  versions  of  this  equation  are 

e 

possible  and  are  referred  to  as  nonlinear  eddy  viscosity  models  or  algebraic  Reynolds  stress  models.  They 
all  have  the  same  fundamental  problem  for  predicting  transition.  These  models  assume  that  the  fluctuations 
are  somehow  in  equilibrium  with  the  mean  flow.  They  suppose  that  a  given  shear  level  results  in  a  given 
turbulence  level.  The  transition  process  is  as  far  from  an  equilibrium  situation  as  possible.  During 
transition,  fluctuation  amplitudes  grow  exponentially  in  time  (moving  with  the  fluid),  while  shear  levels 
remain  almost  constant. 

Despite  these  difficulties,  it  is  interesting  to  note  that  two-equation  models  still  can  produce  qualitatively 
transition-like  behavior,  with  a  rapid  increase  in  turbulence  levels  and  skin  friction  from  low  initial  levels. 
This  suggests  that  these  differential  equations,  and  by  inference  the  modeled  RST  equations,  have  the 
mathematical  capability  of  produce  exponentially  growing  solutions,  and  transitional  behavior.  With  the 
added  accuracy  in  the  production  and  rapid  pressure-strain  terms,  it  is  expected  that  RST  models  could 
potentially  constitute  a  very  accurate  transition  prediction  methodology.  However,  classic  RST  models 
represent  a  significant  investment  in  programming  complexity  and  computational  resources.  The  turbulent 
potential  model  used  in  this  work  (Perot,  1999)  (Perot  &  Wang,  1999)  (Zhang  &  Perot,  2000)  represents  an 
effective  compromise.  It  is  a  reformulation  of  the  RST  equations  that  retains  the  non-equilibrium  and 
energy  redistribution  physics,  but  which  can  be  implemented  at  a  computational  cost  and  programming 
complexity  comparable  to  the  very  popular  two  equation  models.  The  turbulent  potential  model’s 
formulation  and  success  in  predicting  fully  turbulent  flows  has  been  discussed  in  prior  publications  (Perot 
&  Taupier,  1999)  (Tsuei  &  Perot,  2000)  (Are,  Zhang  &  Perot,  2001).  In  this  work  we  focus  on  its  unique 
abilities  to  predict  transition  in  boundary  layer  flows.  We  will  demonstrate  the  ability  to  predict  by-pass 
transition,  the  affects  of  pressure  gradients,  natural  transition,  relaminarization,  and  the  affect  of  noise 
levels  on  natural  transition. 
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3.4  Flat  plate  boundary  layer  in  zero  pressure  gradient 


The  process  of  transition  is  studied  by  looking  at  the  evolution  of  the  friction  coefficient  along  the 
streamwise  direction.  The  friction  coefficient  is  a  very  sensitive  indicator  of  transition  that  increases 
dramatically  as  transition  occurs.  The  model  predictions  are  compared  to  experimental  data  with  different 
turbulence  intensities.  The  mean  velocity  is  initially  uniform  flow  for  all  cases  and  the  initial  values  of 

velocity  U0 ,  turbulence  Reynolds  number  Rer  =  k2/(v£) ,  turbulence  intensity  level 


for  five  experimental  cases  with  different  turbulence  intensities  are  given  in  Table 


3. 1 .  The  initial  values  of  the  turbulent  kinetic  energy  k  are  determined  using 


k  =  3/2(Tu  *£/0)2  *  For  the 


Tu  =  0.03% ,  Tu  =  1 .25%  and  Tu  =  1 .3%  cases,  the  initial  turbulent  dissipation  rate  £  is  calculated 
from  Rer  using  £  =  £2/(vRer)  and  the  value  of  Rer  is  assumed.  Results  are  not  very  sensitive  to 

reasonable  values  of  Rer .  For  the  Tu  -  3%  (T3A)  and  Tu  -  6%  (T3B)  cases,  the  data  for  k  is 

dk 

available,  and  the  initial  values  of  £  are  determined  from  £  —  ,  where  the  subscript  M°°" 

dx 

represents  parameters  in  the  free  stream.  The  values  of  Rer  for  these  two  cases  listed  in  Table  3.1  are  just 
calculated  from  £ .  The  initial  potentials  (j)  and  If/  are  set  as  2/3  k  and  zero  respectively.  All 
experiments  were  performed  in  air  so  V  =  1.55x10  5  was  used  in  every  case. 


U0(m!s) 

Rer 

Tu 

Source 

24.4 

100 

0.03% 

Schubauer  &  Klebanoff 

22 

250 

1.25% 

Abu  -Ghannam  &  Shaw 

14.42 

250 

1.3% 

Dhawan  &  Narasimha 

5.4 

200 

3.0% 

ERCOFTAC,  T3A 

9.4 

200 

6.0% 

ERCOFTAC,  T3B 

Table  3.1  Initial  flow  parameters  at  the  leading  edge  of  a  zero  pressure  gradient  boundary  layer 

The  friction  coefficient  is  plotted  against,  Rex ,  the  Reynolds  number  based  on  downstream  position,  in 
Figure  14.  The  friction  coefficient  correlations  for  laminar  and  turbulence  flows,  Cj  =  0.664 Re~^2  and 
Cf  =  0.027  Re^7  respectively,  are  also  plotted  with  dashed  lines  for  comparison.  Overall,  the  ability  of 

the  present  model  to  predict  the  critical  transition  location  is  quite  good.  However,  the  natural  transition 
case,  Tu=0.03%,  shows  some  discrepancy  with  the  experimental  data  of  Schubauer  (1955)  at  2.8xl06. 
Actually,  this  discrepancy  is  somewhat  expected  given  the  model  initial  conditions.  Experiments  of  natural 
transition  show  a  wide  range  of  transition  locations  from  a  high  value  of  Rex  =  5.0 xlO6  measured  by 
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Wells  (1967)  to  classical  predictive  theories  such  as  e9  rule  which  predict  a  value  of  2.0  xlO6  (Warsi, 
1999).  This  range  of  results  is  commonly  attributed  to  different  noise  levels  in  the  various  experiments. 
We  show  later  that  the  model  can  actually  predict  the  entire  range  of  natural  transition  experiments  by 
varying  the  model  initial  conditions.  We  believe  that  the  ability  of  this  model  to  predict  natural  transition  is 
a  unique  characteristic  that  is  not  demonstrated  by  other  RANS  models.  For  the  highest  inlet  turbulence 
intensity  case  (T3B,  Tu=6.0%),  the  level  of  friction  coefficient  turns  out  to  be  at  high  values  in  the  pre- 
transitional  region.  A  recent  DNS  study  (Jocobs  and  Durbin,  2001)  also  reported  a  similar  over  prediction 
for  a  Tu=7.0%  case. 


Figure  14.  Transition  in  zero  pressure  gradient  boundary  layer  at  various  initial  turbulence  intensities. 
The  symbols  represent  experiment  data  and  the  lines  are  the  present  results.  Dashed  lines  are  fully 
laminar  and  fully  turbulent  correlations. 

The  other  significant  departure  of  the  model  predictions  from  the  experiments  is  in  the  overshoot  in  the 
friction  coefficient  that  occurs  at  the  very  end  of  the  transition  process  and  in  the  transition  length. 
Predicting  this  overshoot  and  length  correctly  is  not  as  important  as  predicting  the  transition  location  itself. 
Interestingly,  the  model  does  show  the  characteristic  overshoot  behavior  found  in  experiments,  but  over 
predicts  the  extent  of  this  overshoot.  Our  current  hypothesis  is  that  this  discrepancy  is  largely  caused  by 
the  use  of  boundary  layer  equations  to  solve  for  the  velocity  and  turbulence  variables,  and  is  not  due  to  the 
model  itself.  The  boundary  layer  equations  are  based  on  the  premise  that  streamwise  second  derivatives 
are  small.  During  the  later  part  of  transition,  as  the  flow  suddenly  becomes  turbulent,  the  boundary  layer 
grows  quite  dramatically,  and  the  assumption  of  small  streamwise  second  derivatives  is  not  very  accurate. 
Including  these  streamwise  derivatives  is  expected  to  add  considerable  diffusion  and  significantly  reduce 
the  current  excessive  overshoot  behavior.  The  fact  that  the  overshoot  error  increases  as  the  turbulence 
intensity  decreases  is  consistent  with  this  behavior,  because  the  lower  intensity  transition  events  are  also 
more  abrupt  and  therefore  violate  the  boundary  layer  approximations  the  most  severely. 

Two  other  numerical  predictions,  both  performed  by  Suzen  et.  al.  (2000),  but  using  their  k  —  £ 
intermittency  model  and  another  k  —  £  model  (Launder  &  Sharma,  1974)  respectively,  are  plotted  in 
Figure  15  for  the  T3A  case.  The  k  -  £  model  predictions  are  not  particularly  good  for  this  particular  case. 
In  contrast,  Suzen’s  model  captured  the  late  stage  of  the  transition  region  (including  the  overshoot)  very 
well,  but  the  transition  onset  was  delayed.  The  present  calculation,  as  shown  in  Figure  15,  shows  a 
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smoother  transition  behavior.  In  other  words,  the  present  model  predicts  not  only  a  correct  position  of 
transition  onset,  but  also  a  fairly  accurate  length  of  the  transition  region. 


Figure  15.  Comparison  of  predictions  of  transition  in  zero  pressure  gradient  boundary  layer  for  T3A  case. 
Symbols  represent  experiment  data  and  lines  are  the  predictions. 

The  comparison  of  the  Reynolds  number  based  on  momentum  thickness,  Re  e  ,  and  shape  factor,  H,  with 

experiments  and  results  from  other  numerical  models  are  shown  in  Figure  16  and  Figure  17  respectively  for 
T3A  case.  It  can  be  noted  that  the  present  prediction  always  stays  between  the  other  two  simulations.  The 
fact  that  Suzen’s  model  has  a  slight  later  prediction  of  the  onset  of  transition  is  also  visible  from  the  H 
curves.  The  streamwise  mean  velocity  profiles  calculated  using  the  present  model  are  compared  at  four 
downstream  stations  with  the  experimental  T3A  data  and  Suzen's  predictions  in  Figure  18(a)  through 
Figure  18(d).  It  is  evident  from  these  figures  that  the  model  predictions  agree  very  well  with  the 
experimental  data  even  for  the  detailed  mean  velocity  profiles. 
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Figure  17.  Comparison  of  H  for  T3A  case. 

3.5  Flat  plate  boundary  layer  in  non-zero  pressure  gradient 

Model  predictions  of  transition  in  two  variable  pressure  gradient  boundary  layers  are  compared  with  the 
experimental  data  of  Coupland  (1990)  in  Figure  19.  The  initialturbulence  levels  of  these  two  cases  are 
3.0%  (T3C3)  and  6.6%  (T3C1)  respectively.  For  both  of  these  flows,  the  pressure  gradient  is  initially 
negative  (favorable)  and  then  positive  (adverse)  in  a  profile  that  was  designed  to  roughly  approximate  the 
flow  over  a  turbine  blade.  The  pressure  gradient  profiles  for  the  two  different  cases  are  roughly  equivalent. 
The  initial  conditions  are  listed  in  Table  3.2  and  are  fully  determined  by  the  initial  experimental  data.  It  is 
very  difficult  to  use  correlations  to  predict  this  variable  pressure  gradient  flow  because  the  pressure 
gradient  does  not  correspond  closely  to  any  single  Falkner-Skan  situation.  In  addition,  the  RANS  models 
tested  by  Savill  (1990)  showed  some  difficulty  with  these  cases. 

The  model  predicts  the  transition  location  in  these  complex  variable  pressure-gradient  boundary  layers 
quite  well.  Note  that  in  both  cases  the  initially  favorable  pressure  gradient  has  delayed  transition  compared 
to  the  flat  plate  case.  The  pressure  gradient  data  is  not  provided  beyond  the  range  of  the  experimental  data, 
so  we  can  not  continue  the  3%  case  beyond  the  experimental  data.  Predictions  given  by  Suzen  et.  al.  (2000) 
are  compared  with  the  present  calculation  for  T3C1  case  again  in  Figure  20. 


Figure  18.  Mean  streamwise  velocity  profiles  for  T3a  case,  (a)  Re,  =  1.35xl05,  (b)  Re,  =2.04xl05, 
(c)  Re,  =  2.74x10s,  (d)  Re,  =4.19xl05. 
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U0  (m/s) 

Rer 

Tu 

5.9 

160 

6.6% 

ERCOFTAC,  T3C1 

3.7 

100 

3.0% 

ERCOFTAC,  T3C3 

Table  3.2  Initial  flow  parameters  at  the  leading  edge  of  a  variable  pressure  gradient  flat  plate  boundary 
layer 


Rex 


Figure  19.  Transition  in  non-zero  pressure  gradient  boundary  layer  at  various  initial  turbulence  intensities. 
The  symbols  represent  experiment  data  and  the  lines  are  the  present  results. 


0.01 


Rex 


Figure  20.  Comparison  of  predictions  of  transition  in  non-zero  pressure  gradient  boundary  layer  for  T3C1 
case.  Symbols  represent  experiment  data  and  lines  are  the  predictions. 

Figure  21  shows  the  comparison  of  Re^  computed  by  the  present  model  and  the  two  other  models  with  the 

experiment  data  for  T3C1  case.  Predictions  of  the  shape  factor  profiles  using  the  three  models  are  plotted  in 
Figure  22  along  with  the  experimental  data.  Again,  the  present  model  exhibits  excellent  performance  in 
calculating  the  shape  factor  in  the  transition  region. 


Re, 

Figure  21.  Comparison  of  Re  e  for  T3C  lease. 
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Figure  22.  Comparison  of  H  for  T3C1  case. 

3.6  Acoustic  Effects  on  Natural  Transition 

It  was  mentioned  previously  that  acoustic  effects  (wind  tunnel  noise)  can  have  a  fairly  strong  influence  on 
natural  transition.  The  experiment  of  Wells  (1967)  is  commonly  assumed  to  be  largely  free  of  acoustic 
noise,  and  therefore  gives  a  high  "quiet"  transition  location  of  Re^  tr  =5.0xl06.  Our  initial  prediction  of 

Rej/r  =  2.0  xlO6  is  close  to  the  low  (noisy)  end  of  the  experimental  range.  The  experiments  of  Schubauer 

(1955)  attempted  to  reduce  noise  levels,  but  were  clearly  not  as  successful  in  this  respect  as  the  later 
experiments  of  Wells. 

Interestingly,  the  turbulent  scalar  potential  represents  irrotational  fluctuations  (see  Perot  1 999),  and  could 
be  a  very  good  measure  of  the  average  acoustic  forcing.  Our  initial  natural  transition  prediction  (Fig.  2.1) 
was  with  relatively  high  noise  levels  of  0  =  2/3  k,  and  is  actually  quite  close  to  “noisy”  experiments.  If  one 
reduces  the  initial  scalar  potential  0  to  a  value  of  0.12k  ,  the  onset  of  the  transition  will  be  delayed  to 

5.0xl06,  which  agrees  with  Wells  (1967).  One  can  also  match  the  data  of  Schubauer  (1955) 
Rcx  (r  =  2.8xl06  in  Table  2.1  by  using  0  =  0.4A: .  AH  three  simulations  are  shown  in  the  Cf  -Rex  curve 

plotted  in  Figure  23.  This  demonstrates  the  models  ability  to  accurately  capture  the  affects  of  noise  on 
natural  transition.  The  effects  of  noise  on  bypass  transition  are  present  but  are  less  significant. 
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Figure  23.  The  acoustic  effects  on  the  natural  transition  at  Tu  =0.03%.  Dashed  lines  are  fully  laminar 
and  fully  turbulent  correlations. 

3.7  Relaminarization  in  Turbulent  Channel  Flow 

Almost  as  critical  as  transition,  is  the  process  of  relaminarization,  where  a  turbulent  flow  decays  to  a 
laminar  state.  Relaminarization  can  occur  because  of  strong  stratification,  strong  rotation,  or  in  the  case  of 
channel  flow,  due  to  a  reduction  in  the  driving  pressure  gradient.  Figure  24  shows  predictions  of  the  skin 

friction  in  channel  flow  for  various  bulk  Reynolds  numbers  Re  =  ,  where  h  is  the  channel  half  height 

V 

and  Ub  is  the  average  velocity.  The  figure  also  has  dashed  lines  for  the  laminar  exact  result  Cj  =  |Re~l 
and  the  turbulent  correlation  of  Dean  (1978)  Cf  =  0.044  Re~°  227 .  Here  the  skin  friction  coefficient  is 

T 

calculated  using  Cf  = - liL~2y  where  Tw  is  the  wall  shear  stress  and  u0  is  the  center  line  velocity  of  the 

lPU 0 

channel.  The  model  shows  a  clear  relaminarization  at  Re6  =  1700  (Rer  =  120).  Experiments  and  DNS 
(Keefe,  1992)  indicate  a  transition  Reynolds  number  of  close  to  1000  (Rer  =  60 ).  So  the  model  could  be 
improved.  However,  a  theoretical  study  by  Orszag  (1971)  showed  the  linear  instability  for  a  channel  flow 
occurs  at  Re^,  =  1600  (Rer  =110).  It  is  our  observation,  that  this  relaminarization  location  is  strongly 

influenced  by  the  low  Reynolds  number  behavior  of  the  model  (particularly  the  dissipation  and  pressure 
strain  term),  and  not  the  transition  term  itself. 
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Figure  24.  predictions  of  the  skin  friction  in  channel  flow  for  various  bulk  Reynolds  numbers.  Dashed  lines 
are  fully  laminar  and  fully  turbulent  correlations. 


3.8  Conclusions 

The  ability  of  the  turbulent  potential  model  to  predict  transition  in  a  wide  variety  of  boundary  layer  flows 
has  been  demonstrated.  This  includes  the  ability  to  predict  everything  from  natural  transition  (a  first  for 
RANS  based  models)  to  large  free-stream  turbulence  intensities.  The  model  is  also  demonstrated  to  be  able 
to  accurately  calculate  the  transitional  flat  plate  boundary  layer  under  favorable  and  adverse  pressure- 
gradients.  The  ability  of  this  model  to  reflect  the  effect  of  noise  on  the  natural  transition  is  also  shown  by 
comparison  with  corresponding  experimental  data  under  different  conditions.  Finally,  the  relaminarization 
in  turbulent  channel  flow  is  successfully  captured  using  this  model.  The  computation  cost  of  the  proposed 
method  is  comparable  to  the  widely  used  two-equation  models  such  as  k  /  £ .  However,  unlike  k  /  £ ,  the 
model  does  not  assume  equilibrium  between  the  developing  instability  and/or  turbulence  fluctuations,  and 
it  is  this  non-equilibrium  nature  to  which  we  attribute  these  successful  predictions. 


4.  Large  Eddy  Simulation 

This  portion  of  the  research  demonstrates  that  the  turbulent  potential  model  can  be  used  to  model  small 
scale  turbulence  in  the  context  of  a  large  eddy  simulation.  This  is  essentially  using  a  RANS  model  to  do 
LES.  But  in  this  context  we  must  be  very  careful  to  delineate  what  is  meant  by  a  RANS  model  and  LES 
because  the  formal  definitions  of  these  two  modeling  approaches  makes  them  mutually  exclusive.  What 
we  intend  by  the  statement  ‘RANS  model  in  an  LES  context’  is  the  use  of  models  which  involve  the 
solution  of  partial  differential  equations  to  obtain  the  unknowns  necessary  to  close  and  solve  the  system. 
Standard  LES  subgrid  scale  models  are  algebraic  in  nature. 

The  use  of  RANS  models  in  an  LES  context  is  not  immediately  intuitive  from  either  the  RANS  or  LES 
perspective.  LES  subgrid  models  have  always  had  a  focus  on  simplicity  and  low  cost,  and  the  solution  of 
partial  differential  equations  is  rarely  is  considered  to  be  either.  The  construction  of  classic  RANS  models 
relies  heavily  on  notions  of  turbulence  universality  and  often  turbulence  equilibrium,  both  of  which  are 
harder  to  justify  if  we  now  restrict  the  model  to  represent  only  fluctuations  below  a  subgrid  threshold. 
Nevertheless,  we  intend  in  this  section,  to  show  the  benefits  of  using  just  this  counterintuitive  approach.  It 
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will  be  shown  that  some  RANS  models  (such  as  the  turbulent  potential  model)  can  be  applied  successfully 
well  beyond  their  formal  range  of  applicability  and  into  the  subgrid  modeling  regime. 

The  use  of  RANS  models  to  perform  large  eddy  simulation  is  not  an  entirely  new  idea.  Deardorf  (1973) 
used  a  Reynolds  stress  transport  equation  model  in  a  LES  of  atmospheric  turbulence.  Soon  after  this, 
Schumann  (1975)  performed  a  large  eddy  simulation  using  a  partial  differential  equation  for  the  unresolved 
turbulent  kinetic  energy.  More  recently  Glosal,  Lund,  &  Moin  (1995)  have  tested  a  one-equation  LES 
model  which  also  uses  a  kinetic  energy  evolution  equation.  They  concluded  that  the  model  worked  well, 
but  that  the  additional  cost  incurred  was  not  advantageous  compared  to  traditional  LES  models.  Fureby  et 
al  (1997)  have  essentially  reproduced  Deardorf  s  earlier  work  using  updated  closure  models,  and  have 
shown  excellent  results  for  isotropic  decay  and  turbulent  channel  flow  with  only  a  20%  increase  in 
computational  cost  over  standard  LES  models.  Speziale  (1998)  proposed  a  switch  for  the  k/e  model  that 
would  allow  it  to  smoothly  transition  into  an  LES  subgrid  scale  model.  Along  similar  lines,  Spallart  (2000) 
has  demonstrated  a  switch  for  the  Spallart-Allmaras  (1992)  model  that  allows  it  to  smoothly  transition  from 
RANS  to  LES.  The  hybrid  Spallart  model  is  called  detached  eddy  simulation  (DES)  because  the  switch  is 
based  on  the  distance  to  the  wall,  and  therefore  causes  RANS  type  behavior  near  the  wall  and  LES  far  away 
(for  the  detached  eddies).  A  number  of  simulations  based  on  both  the  Speziale  and  DES  approaches  have 
subsequently  been  performed  (Peltier  et  al.  2000,  Travin  et  al  2000,  Arunajatesan  et  al.  2000). 

The  fundamental  difference  between  the  proposed  approach  and  the  hybrid  models  of  Speziale  and  Spallart 
is  that  the  turbulent  model  does  not  require  an  explicit  switch  to  do  LES.  This  is  actually  a  profound 
difference,  but  the  importance  of  the  distinction  is  very  difficult  to  convey  abstractly. 


4.1  Theoretical  Approach 

We  start  with  the  decomposition  of  the  instantaneous  velocity  field  Vi : 

V,=U,{A,t)+u,  (4.D 

where  £/.  and  uf  are  the  resolved  and  unresolved  (modeled)  velocities  respectively.  The  resolved  part  Ui 
accounts  for  the  flow  motion  of  unsteady,  energy  containing  scales  and  therefore  is  a  function  of 
wavenumber  X  and  time  t .  It  also  accounts  for  any  steady,  mean  motion.  The  evolution  equation  for  the 
resolved  velocity  field  can  be  written  as 

+  =  +  (4.2, 

dt  1  dxj  dxj  dxjdxj 

where  the  angular  bracket  represents  averaging  over  modeled  scales.  The  resolved  pressure  pr  is  defined 
by 


»,  SU,  dUJ  |  (4,3) 

dxjdxi  dxj  dXj  dxjdxj 

The  primary  variable  to  be  numerically  solved  in  this  work  is  the  stream  function,  so  the  pressure  is  only 
explicitly  required  at  boundary  cells.  The  sum  of  the  last  term  on  the  r.h.s.  of  Equation  (4.2)  can  be 
considered  playing  the  same  role  as  the  convection  term  in  the  filtered  momentum  equation  of  a 
conventional  LES  method.  Instead  of  using  the  filtering  function  to  account  for  the  contribution  of  the 
Reynolds  stress  tensor  (Residual  stress  tensor)  in  the  momentum  equation,  we  stick  to  relating  it  to  the 
turbulent  potentials.  So  the  equation  (3.2)  is  rewritten  as: 
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(4.4) 


Wl{UWt_  »r||.  dWk 

dt  J  dxj  9x,  dxjdxj  >jk  dxj  9x,. 

The  flow  motion  above  the  grid  scale  should  be  resolved  by  the  numerical  solver,  and  evolve  via  equation 
(4.4),  which  can  be  closed  by  the  TPM  model  solving  for  the  small  scale  motions.  The  model  equations  can 
be  found  in  the  work  of  Zhang  (2002)  with  a  few  modifications. 

The  transport  equation  for  the  resolved  kinetic  energy  can  be  obtained  by  taking  the  dot  product  of  (4.4) 
with  the  velocity  vector. 


dKr  TT  dKr  d  trr  \  3 

~ir+u’^=-^{u‘p')+y^ 


u , 


W L 

'  dxJ  J 


OX; 


-v- 


dUj  dU, 

dxj  dxj 


-£:n 


Uk 


dxj 


kU; 


(4.5) 


The  evolution  of  the  total  resolved  kinetic  energy  can  be  derived  by  the  integrating  equation  (4.5)  over  the 
computational  domain.  Notice  that  the  convection  term  on  the  l.h.s.  and  the  first  three  terms  on  the  r.h.s. 
just  redistribute  the  kinetic  energy  inside  the  periodic  domain.  We  end  up  with  a  compact  form: 


J^(j Krdv)=-\£rdv-  \Prdv 


(4.6) 


7)1  J  7)1/ 

where  the  resolved  dissipation  term  is  e  -y  .J . — L  and  the  resolved  production  term  is 


dxj  3 Xj 


P.=eM 


y* 


*Yk 

dxj 


Ui .  While  this  is  called  a  production  term  because  of  its  form,  this  term  tends  to  reduce 


the  resolved  kinetic  energy  (sending  it  to  the  unresolved  kinetic  energy).  The  evolution  equation  for  the 
modeled  kinetic  energy  was  given  by  (4.5).  The  total  modeled  kinetic  energy  is  governed  by  the  following 
integral  equation, 


^(f  KJV)=-\£JV+  \PmdV 


(4.7) 


where  the  resolved  production  term  is  Pm  =  (O  Xf/l  and  Q)i  is  a  component  of  vorticity  vector.  Note  that 

the  production  terms  for  the  resolved  and  modeled  energy  equations  are  mathematically  equivalent,  though 
they  were  calculated  in  the  CFD  code  in  the  different  forms.  This  can  be  more  easily  seen  from  simple 
derivation, 


\prdv=  \eijkd^uidv=  l^-(£ijk¥kUj)dv- Kfr-iMv 


dxj 


=  Ykdv=  \«>k¥kdv 


(4.8) 
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4.2  Channel  Flow 


Shown  below  are  various  calculations  using  the  turbulent  potential  model  of  turbulent  channel  flow  at  a 
shear  velocity  Reynolds  number  of  1 80.  The  channel  flow  is  a  very  well  understood  flow,  both  in  terms  of 
LES  models  and  DNS  statistics  of  the  flow.  The  low  Reynolds  number  was  chosen  so  that  we  could  test 
the  model  right  into  the  DNS  limit.  The  domain  is  periodic  in  the  streamwise  and  spanwise  directions  and 
only  half  the  wall  normal  domain  is  simulated  with  a  symmetry  condition  at  the  channel  midplane.  The 
result  is  therefore  closer  to  an  open  channel  flow  simulation,  but  it  is  observed  (Salvetti  et  al ,  1997)  that 
this  midplane  boundary  condition  does  not  affect  the  wall  boundary  layer  turbulent  statistics. 

Five  cases  were  run  with  varying  mesh  resolutions  in  the  spanwise  and  streamwise  directions.  The  wall 
normal  resolution  of  60  points  is  required  even  for  a  pure  RANS  simulation  of  the  channel  if  wall  functions 
are  not  used,  so  the  wall  normal  resolution  was  left  unaltered.  The  simulations  were  initialized  with 
random  fluctuations  and  a  uniform  mean  velocity  in  the  streamwise  direction.  The  simulation  was  then  run 
until  statistical  steady  state  was  achieved.  A  series  of  five  simulations  were  run  with  all  variables  fixed 
except  the  spanwise  and  streamwise  mesh  resolution.  Table  3.1,  shows  the  resolutions  for  the  various 
simulations  and  the  resolution  of  the  reference  DNS  simulation  of  Moser,  Kim  &  Mansour  (1999).  The 
current  simulations  vary  from  DNS  resolution  to  LES  type  resolutions  to  what  should  closely  approximate 
RANS  resolutions  (where  the  near  wall  streaks  can  not  possibly  be  resolved).  Case  1  used  an  80x60x80 
Cartesian  mesh.  Each  succeeding  case  doubles  the  streamwise  and  spanwise  grid  size  so  that  Case  5  uses  a 
5x60x5  grid.  Figure  25  shows  the  mean  velocity  profiles  for  all  the  simulations  and  the  DNS  results  of 
Moser  et  al. 


Moser,  Kim 
&  Mansour 

Case  1 

Case  2 

Case  3 

Case  4 

Case  5 

Ax+ 

17.7 

18.0 

36.0 

72.0 

144.0 

288.0 

EM 

0.054 

0.48 

0.48 

0.48 

0.48 

0.48 

A z+ 

5.9 

6.75 

13.5 

17.0 

34.0 

68.0 

Ax//? 

0.098 

0.10 

0.20 

0.40 

0.80 

1.60 

(4y/^)m  ax 

0.024 

0.014 

0.014 

0.014 

0.014 

0.014 

Az/h 

0.033 

0.037 

0.075 

0.15 

0.30 

0.60 

Table  3.1:  Mesh  resolutions  for  calculations  of  channel  flow  using  the  turbulent  potential  model. 

The  important  result  demonstrated  in  figure  25  is  that  the  turbulent  potential  model  naturally  adapts  to  any 
mesh  resolution.  The  skin  friction  is  within  3%  of  the  DNS  data  for  all  cases.  Case  1  is  essentially  DNS, 
the  modeled  turbulence  quantities  are  all  extremely  small  in  this  case.  Case  2  is  LES,  the  turbulent  eddy 
viscosity  is  roughly  the  order  of  the  laminar  viscosity  in  this  case.  (Note  that  eddy  viscosity  is  not  used  in 
the  turbulent  potential  model  but  was  calculated  for  these  test  cases  as  a  simple  way  of  evaluating  the 
importance  of  the  model  on  the  resolved  flow).  Case  3  shows  the  most  model  overlap,  the  flow  solution  is 
still  unsteady,  but  the  flow  is  probably  dominated  by  the  model  at  this  point.  The  eddy  viscosity  is  up  to 
ten  times  the  laminar  viscosity,  and  the  resolved  RMS  levels  are  only  a  third  of  their  DNS  levels.  Cases  4 
and  5  are  fully  RANS  simulations  and  display  no  unsteadiness  or  three-dimensionality  at  steady  state. 
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Figure  25.  Calculations  of  mean  velocity  in  a  turbulent  channel  flow  atRer  =180  using  the  turbulent 
potential  model.  The  mesh  resolutions  ranging  from  DNS  to  RANS,  and  are  compared  with  the  DNS  data 
of  Moser  et  al . 


Note  that  hybrid  models  that  are  based  on  a  switch  (al!  existing  ones)  can  not  reproduce  Figure  25.  The 
explicit  switch  invariably  alters  the  model  in  the  coarse  mesh  limit  resulting  in  mesh  dependent  results  (see 
Nikitin  et  al.  2000).  It  is  believed  that  the  result  presented  above  is  closely  tied  to  the  non-equilibrium 
nature  of  the  model.  Tests  of  a  low  Reynolds  number  k/e  model  (not  shown)  indicated  that  k/e  could  not 
function  in  the  LES  and  DNS  regime  unless  the  explicit  Speziale  damping  switch  was  imposed.  This 
switch  invariably  altered  the  coarse  mesh  solutions.  It  is  hypothesized  that  Reynolds  stress  transport 
equations  could  reproduce  Figure  25,  but  given  the  reduced  cost  and  increased  robustness  of  the  turbulent 
potential  model,  the  advantages  of  using  a  Reynolds  stress  transport  model  are  limited. 

A  preliminary  observation  that  was  obtained  from  these  simulations  is  that  the  transition  from  LES-like 
(large  scale  unsteadiness)  to  RANS-like  (steady,  one  dimensional)  solutions  was  relatively  abrupt.  This 
may  be  related  to  the  low  Reynolds  number  of  the  simulations,  but  Nikitin  et  al.  report  similar  behavior 
with  DES  at  much  larger  Reynolds  numbers  and  suggests  that  the  phenomena  occurs  when  the  mesh  no 
longer  can  resolve  the  larger  scale  energy  containing  turbulent  scales. 

The  turbulent  potential  model  (like  other  hybrid  approaches)  solves  the  issue  of  how  to  perform  LES  near 
walls  at  high  Reynolds  numbers.  Standard  LES  requires  that  turbulence  producing  structures  be  resolved. 
Near  the  wall  the  important  structures  are  streaks  which  scale  inversely  with  the  Reynolds  number.  This 
means  that  standard  LES  of  near  wall  turbulence  which  ought  to  be  Reynolds  number  independent,  is  not. 
Significant  efforts  to  generate  wall  functions,  or  discontinuous  two  layer  models  to  avoid  this  Reynolds 
number  dependence  have  lead  to  methods  for  the  wall  which  display  a  great  deal  of  complexity  (Cabot  & 
Moin  1999,  Baggett  1997,  Bagwell  et  al.  1993).  LES  models  which  have  the  ability  to  transition  smoothly 
to  a  RANS  limit  (such  as  the  proposed  approach)  can  model  near  wall  turbulence  with  only  a  possible 
logarithmic  dependence  on  the  Reynolds  number. 

Finally  it  should  be  noted  that  the  fact  that  the  proposed  model  can  be  applied  in  both  the  RANS  and  LES 
limits  makes  it  extremely  attractive  to  the  end  users  of  such  models.  The  ability  to  use  a  whole  range  of 
modeling  solutions  within  a  single  code  and  unified  modeling  framework  greatly  enhances  the  users  ability 
to  make  the  difficult  choice  between  cost  and  accuracy  that  turbulence  invariably  imposes. 
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4.3  RANS  Modeling  of  Isotropic  Decaying  Turbulence 

All  non-equilibrium  models  should  work  for  isotropic  turbulence  before  it  is  applied  to  other  flows.  A 
uniform  643  quadrilateral  mesh  is  placed  on  a  18/r3  cubic  domain.  The  domain  is  periodic  in  the 
streamwise,  spanwise  and  lateral  directions.  We  first  tested  the  performance  of  the  model  in  predicting  the 
evolution  of  turbulence  quantities  during  the  decay.  The  initial  condition  is  taken  from  a  5 12 3  simulation 
previously  calculated  by  de  Bruyn  Kops  (1998).  This  case  is  equivalent  to  the  famous  decay  experiment  of 

Comte-Bellot  and  Corrsin  (1971).  The  flow  has  an  initial  turbulent  kinetic  energy  of  753  (cm2  / S2  )  and 

the  initial  dissipation  rate  of  the  value  4487  ( cm2 / s 3  ).  This  initial  condition  matches  the  laboratoiy 
kinetic  energy  at  xl M  =  42  where  M  is  the  grid  space  in  the  experiment. 

The  model  was  tested  under  the  condition  that  the  resolved  velocity  is  negligible.  This  is  a  full  RANS 
simulation,  with  no  perceptible  resolved  flow,  even  though  the  mesh  is  643 .  The  prediction  of  present 
model  was  shown  in  Figure  26  and  it  is  consistent  with  the  DNS  data  of  de  Bruyn  Kops. 


Figure  26.  RANS  simulation 


4.4  LES  of  Isotropic  decaying  turbulence 

This  simulation  took  a  643  truncation  of  the  5123  velocity  field  computed  by  de  Bruyn  Kops  (1998)  as  the 
initial  condition  for  the  resolved  flow  field.  This  velocity  field  was  generated  in  a  series  of  process  so  that 
the  initial  energy  spectrum  matches  the  experiment  and  the  large-scale  structures  are  properly  developed. 
The  contour  plot  of  the  initial  resolved  velocity  component  is  shown  in  Figure  27. 
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Figure  27.  Velocity  contour  of  initial  isotropic  turbulence 

It  can  be  calculated  that  the  initial  values  for  the  resolved  kinetic  energy  and  initial  resolved  dissipation  are 
552  ( cm2 / S 2  )  and  395  ( cm2 / S 3  )  respectively.  To  match  the  total  initial  kinetic  energy  and  initial 
dissipation  rate  of  the  experiment,  the  initial  values  of  the  modeled  turbulence  quantities  are  taken  to  be 
201  {cm2/ S2  )  for  the  kinetic  energy  and  4192  {cm2 f  s*  )  for  the  dissipation  rate.  The  decay  of  the 
resolved,  modeled  and  total  kinetic  energies  are  plot  in  Figure  28  together  with  the  result  from  de  Bruyn 
Kops  (1998)  for  comparison.  The  variation  of  time  scale  k! €  is  plotted  in  Figure  29.  For  all  cases,  the 
time  variable  was  non-dimensionalized  with  respect  to  the  initial  eddy  turnover  time  T0  =  K0  /  £0 . 


>, it . 

Figure  28.  Comparison  of  the  kinetic  energy  decaying 


36 


The  decay  of  the  total  kinetic  energy  agrees  well  with  the  DNS  data  in  the  early  and  end  stages  of  the 
simulation.  The  reason  for  the  discrepancy  in  the  middle  range  is  under  investigation.  It  was  found  that  the 
decay  rate  of  the  resolved  scales  is  affected  significantly  by  the  way  it  communicates  with  the  modeled 
scales  in  the  actual  CFD  code.  So  hopefully  improvement  in  the  numerical  scheme  for  the  resolved  flow 
field  will  produce  a  more  reasonable  solution. 


Figure  30.  The  fracture  of  kinetic  energies 


Figure  29  should  be  linear.  The  results  of  DNS  and  the  RANS  model  are  indeed  linear.  The  LES  simulation 
is  not  linear  and  is  the  focus  of  the  investigation.  Figure  30  is  the  ratios  of  modeled  and  resolved  kinetic 
energies  to  the  total  kinetic  energy.  It  can  be  observed  that  the  ratio  of  the  modeled  portion  is  increasing 
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rapidly.  This  implies  that  the  modeled  scales  might  be  excessively  energetic.  The  similar  trend  was  reported 
by  Girimaji  (2002)  and  considered  unphysical. 

Figure  3 1  and  Figure  32  show  the  evolution  of  production  and  dissipations  respectively.  The  resolved  and 
modeled  production  terms  match  each  other  exactly  as  expected.  The  magnitude  of  the  modeled  dissipation 
is  much  bigger  than  its  resolved  counterpart  but  the  latter  wins  quickly  over  the  later  stage.  The  resolved 
dissipation  also  starts  to  decrease  in  magnitude  soon  after  that. 


Figure  32  Evolution  of  Dissipation 
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The  derivatives  of  resolved  and  modeled  kinetic  energies  are  plotted  in  Figure  33  with  the  evolution  of 
-J £rdv~  J Prdv  and  ~^£mdv  +  jPmdv  (the  rhs  of  eQuation  4-6  and  4  7  respectively).  From  this 

picture,  it  can  be  seen  that  the  conservation  of  energy  was  correctly  captured  by  the  numerical  computation. 
It  is  not  clear  whether  the  resolved  dissipation  was  calculated  accurately  enough. 


Term 

magnitude 


Figure  33.  Conservation  of  energy  equation  terms 


4.5  Conclusion 

The  simulation  results  have  shown  that  the  turbulence  potential  model,  when  implemented  as  a  pure  RANS 
model,  is  able  to  produce  predictions  in  accordance  with  DNS  for  the  isotropic  decaying  turbulence.  When 
the  model  is  used  as  a  LES  model,  the  results  are  fairly  close  to  the  DNS  data,  but  could  use  further 
improvement. 

The  modeled  scales  tend  to  be  over  energetic  in  the  late  stage.  The  reason  is  under  investigation.  A  control 
parameter  relating  to  the  proportion  of  the  modeled  kinetic  energy  is  considered  to  be  used  to  yield  the 
desired  solution 

An  important  issue  in  LES  is  the  correct  transfer  of  kinetic  energy  between  the  resolved  and  modeled 
(residual)  flow  field.  Globally,  energy  is  transferred  from  the  large  scales  to  the  small  scales.  Nevertheless, 
the  production  in  the  modeled  scales  can  be  negative  locally  and  backscatter  transfer  can  happen.  The 
modeled  initial  kinetic  energy  will  be  set  a  large  value  while  the  resolved  part  is  initially  small.  We  will 
expect  the  resolved  energy  can  be  fed  by  the  decaying  modeled  part  and  rise  up  to  a  certain  value. 

5.  Conclusion 

This  three  year  project  to  evaluate  the  turbulent  potential  model’s  ability  to  predict  rotation,  turbulent 
transition,  and  sub-grid  scale  turbulence  has  been  very  successful.  The  correct  mathematical  framework  for 
implementing  rotation  was  determined,  a  wide  variety  of  boundary  layer  transition  cases  were  predicted, 
and  the  possibility  of  using  the  model  for  LES  was  explored.  In  summary,  it  has  not  only  been 
demonstrated  that  the  turbulent  potential  model  is  competitive  with  existing  RANS  models,  it  has  been 
shown  that  the  model  can  be  applied  to  situations  usually  considered  to  be  outside  the  realm  of  RANS 
models. 
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Appendix 


A1.0  Frame  invariant  quantities  and  frame  consistent  equations 

Frame  consistency  in  the  context  of  turbulence  modeling  can  be  defined  as  the  property  of  the  model  to 
yield  the  same  solution  regardless  of  the  non-inertial  frame  from  which  it  is  being  computed.  Frame 
invariant  quantities  are  quantities  which  do  not  change  from  one  non-inertial  frame  to  another.  It  can  be 
said  that  frame  consistency  is  defined  in  the  context  of  equations,  while  frame  invariance  is  defined  in  the 
context  of  quantities  like  tensor  fields  and  vector  fields. 

We  consider  a  reference  frame  (*)  which  is  in  a  state  of  arbitrary  time  dependent  rotation  with  respect  to  a 
fixed  inertial  frame  and  with  the  same  origin.  The  position  of  a  point  with  respect  to  the  rotating  and  fixed 
frames  are  given  by  x*  and  x  respectively.  Now 

(o.i) 

where  emand  ek  are  the  unit  vectors  representing  the  axes  of  the  fixed  and  the  moving  reference  frame 
respectively.  The  relationship  between  em  and  ek  can  be  given  by 


e>Tkmem  (0.1) 

where  T  is  a  time  dependent  orthogonal  (or  unitary)  transformation  matrix  describing  the  rotation.  Thus, 
from  (0.1)  and  Error!  Reference  source  not  found,  we  can  say  that 

x*=Tx  (0.2) 

If  we  differentiate  this  with  respect  to  time  then  we  get 


v*=Tv+Tx  (0.3) 

Taking  an  average  we  get 

v*=Tv+Tx  (0.4) 

Subtracting  (0.4)  from  (0.3),  and  then  invoking  (0.1),  we  get 

U*=Tu  (0.5) 

which  is  the  definition  of  a  frame  invariant  vector.  From  (0.3)  and  (0.4)  we  can  see  that  the  total  and  mean 
velocity  fields  are  not  frame  invariant.  It  may  be  easier  to  understand  the  meaning  of  frame  invariance  if  we 

take  T=I,  the  identity  matrix,  at  a  certain  instance  of  time.  In  that  case,  V*=V+T  X  and  U*— U  .  Which 
means  that  the  fluctuating  velocity  is  the  same  when  it  is  seen  from  two  frames  overlapped  on  each  other, 
but  the  velocity  transformation  is  dependent  on  the  frame  rotation  rate.  Thus,  to  get  the  components  of  a 
frame  invariant  vector  in  a  different  frame,  we  simply  project  the  vector  onto  axes  of  the  different  frame. 
However  a  vector  which  is  not  frame  invariant  changes  from  one  frame  to  another.  The  Reynolds  stress 
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tensor  can  be  shown  to  be  frame  invariant  (Speziale  (1970)),  though  the  criterion  for  frame  invariance 
changes  since  it  is  a  second  order  tensor: 


R*=TRTt  (0.6) 

The  Reynolds  stress  tensor  transport  equation  has  the  material  derivative  of  the  Reynolds  stress  on  the  left 
hand  side  of  the  equation.  The  material  derivative  of  the  Reynolds  stress  is  not  a  frame  invariant  quantity. 
We  can  show  this  through  the  following  steps,  which  are  written  in  Cartesian  tensor  notation.  From  (0.2) 

x*  =  Tikxk  (0.7) 


Multiplying  both  sides  by  Tim  we  get 


TimXi  =  TimTikXk 


Now  since  T  is  an  orthogonal  matrix,  therefore  TimTik=5rak,  and  hence 


xm=Timxi 


Therefore  we  finally  have  the  relation 


and  the  chain  rule 


T- 


(0.8) 


(0.9) 


(0.10) 


(0.11) 


Similarly,  for  the  time  derivative  we  can  invoke  (0.9)  and  the  fact  that 


t* 


=  t  and  get 


d  _  dt  d  dxk  d 

at*  at*  at  at*  3xk 

_  d  dT^C  d 

dt  3t*  9xk 


=  -  +  T-x„ 


Now,  we  find  out  the  transformation  for  the  material  derivative  operator 


a  a 
at* +  Vk  ax* 


a  •  .  a  .  a 

-^7  +  Tml  Xm^  +TklVk'T~ 

at  dx,  dx. 


(0.12) 


(0.13) 


Using  (0.13)  and  (0.4) 
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(0.14) 


_a_  _♦  a  a  _  _a_ 

3t*  k  3xk  3t  k  3xk 


It  can  be  immediately  seen  that  a  material  derivative  operated  on  a  scalar  quantity  will  be  frame  invariant. 
Now  we  operate  it  on  the  Reynolds  stress  tensor.  Using  (0.14)  and  (0.6)  we  can  get 


9R:  9r:  _  dTimR,Tjn  3R, 

at*  k  ax;  at  im  jn  k  aXk 


=  {Tim  Tj„+Tjm  Tjn  IR^+T^  Rmn  +TimTjnvk 


-  3R„ 


(0.15) 


ax. 


Clearly,  the  material  derivative  is  not  frame  invariant  when  it  is  operated  on  a  second  order  tensor. 

Oldroyd  (1950)  proposed  a  derivative  in  order  to  get  a  time  rate  of  change  of  a  tensor  with  respect  to  a 
coordinate  system  getting  convected  and  distorted  along  with  the  deforming  fluid.  This  gives  a  derivative 
which  is  frame  invariant.  The  Oldroyd  derivative  of  the  Reynolds  stress  tensor  is  given  by 


st  at  kaxk  ikaxk 


-R 


jk 


dvi 

3x, 


(0.16) 


It  can  be  noted  that  the  third  and  the  fourth  terms  are  new  when  compared  to  the  material  derivative.  Using 
steps  similar  to  those  used  from  (0.7)  to  (0.15),  we  can  verify  (Appendix  A)  that 


M=x  t  JK* 


St* 


St 


(0.17) 


where  —  is  being  used  to  denote  the  Oldroyd  derivative.  A  very  important  point  to  be  noted  is  that  the 
8t 

operation  of  an  Oldroyd  derivative  on  a  tensor  gives  a  frame  invariant  quantity  only  if  the  tensor  itself  is 
frame  invariant.  Similar  definition  holds  for  the  Oldroyd  derivative  of  a  vector. 

A1.1  Some  frame  invariant  velocity-field  dependent  quantities 

The  following  velocity-field  dependent  tensors  are  frame  invariant  (Fredrickson  (1964)).  These  are  the 
strain  and  the  intrinsic  vorticity  .  The  strain  tensor  S  is  defined  as 

S  =  I(y  +  v  )  (0.18) 

umn  £  V  m>n  n,m/  v  7 

While  the  intrinsic  vorticity  tensor  W  is  defined  as 

Wmn=Wmn+enmknk  (0.19) 


Here 
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(0.20) 


W  =— (v  —v  ) 

mn  «  \  m,n  n,m ) 


is  the  conventional  vorticity  tensor. 


A1.2  Frame  consistent  equations  for  Reynolds  stresses  and  turbulent  potentials 

If  we  use  an  Oldroyd  derivative  on  the  left  hand  side  of  the  exact  RST  equation,  then  we  get 

Rjj.l  Vk^-ij,k  ^i)  —  ^ij  _  ^ij  ^"ijk.k  ^■^'ij.kk 

Here  we  have  a  new  production  term,  which  includes  the  Coriolis  term 

Py  =  -RikVjk  —  RjkVi>k  —  ejmk^ik^m 
=  -Ru<(Wjk+Sjk)-Rjk(Wik+Sik) 

The  new  production  term  is  a  frame  invariant  quantity  because  it  involves  the  intrinsic  vorticity  tensor  and 
strain  tensor.  On  the  contrary  Py  is  not  a  frame  invariant  quantity.  As  shown  by  Speziale  (1979),  the 
Reynolds  stress,  pressure  strain,  dissipation,  turbulent  diffusion  and  viscous  diffusion  can  also  be  found  to 
be  frame  invariant  quantities. 


(0.21) 


(0.22) 


Thus  the  left  hand  side  of  (0.21)  is  frame  invariant  because  it  is  equal  to  the  Oldroyd  derivative  while  all 
the  terms  on  the  right  hand  side  of  (0.21)  are  frame  invariant,  from  the  discussion  in  the  previous 
paragraph.  There  are  two  things  to  be  inferred  from  this.  Firstly,  the  exact  RST  equation  is  frame 
consistent.  Secondly,  any  Reynolds  stress  transport  model  will  be  material  frame  consistent  if  an  Oldroyd 
derivative  is  used  on  the  left  hand  side,  and  frame  invariant  terms  are  used  on  the  right.  Now  we  go  on  to 
try  and  find  out  a  similar  criterion  for  the  transport  equations  of  the  turbulent  potentials  to  be  frame 
consistent. 

We  first  assume  the  turbulent  potentials  to  be  frame  invariant  quantities.  Then  we  check  whether  this 
assumption  gives  equations  which  have  the  correct  form.  If  that  is  true,  then  turbulence  potentials  are  frame 
invariant  quantities. 

The  assumption  that  turbulent  potentials  are  frame  invariant  means  that 


4>*  =  4> 

Vi  =Tik\|/k 


(0.23) 


and 


av  ■  a  Au 

3x‘kax;  ax’m  ax^ 

=  •  d  cdRL 

dxkdxk  dxq  dxn 


(0.24) 


(0.25) 
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If  the  chain  rules  (0.9)  to  (0.1 1)  along  with  (0.6)  and  (0.23)  are  inserted  in  (0.24)  and  (0.25),  then  we  get 
back  equations  which  have  exactly  the  same  form  as  (0.24)  and  (0.25),  but  with  quantities  which  are  in  the 
original  untransformed  frame  of  reference. 

Since  the  turbulent  potentials  are  frame  invariant,  therefore  the  Oldroyd  derivative  of  the  turbulence 
potentials  will  be  a  frame  invariant  quantity.  Thus  a  reliable  way  to  construct  a  frame  consistent  model 
using  turbulence  potentials  is  to  use  the  Oldroyd  derivative  on  the  left  hand  side  of  the  equation  (transport 
terms),  and  ensure  that  all  the  terms  on  the  right  hand  side  (model  terms)  are  frame  invariant.  The  equation 
for  the  turbulence  potentials  now  take  the  form 

^  +  v-V<])  =  2P$  +  ^  +  V-C$-e,+-uV2(|)  (0.26) 

at 


—  +  v  •  V\|/  -  \|f  Vv  =  2P  +  n  +  V  •  Cv  -  ey  +  vV2\j/  (0.27) 

9t 

Here  the  main  symbols  on  the  right  hand  side  in  (0.26)  and  (0.27)  have  the  same  meaning  as  the  symbols  in 
(0.21),  i.e.  P  denotes  frame  invariant  production,  it  denotes  pressure  strain,  C  denotes  turbulent  diffusion, 
and  £  denotes  the  dissipation.  The  equation  for  <j)  is  a  scalar  equation,  while  the  equation  for  is  a  vector 
equation.  The  Oldroyd  derivative  for  <|>  is  the  familiar  material  derivative,  while  that  for  \|/  is  equal  to 
material  derivative  in  2D  flows  and  for  flows  which  are  inhomogenous  in  only  direction. 

A  small  point  is  to  be  noted  when  the  Oldroyd  derivative  of  R12  and  y3  is  compared.  The  Oldroyd 
derivative  of  R]2  is 


5R12  _  9R12  9R,2  „  9v, 

5t  9t  2  dx2  22  dx2 

While  the  Oldroyd  derivative  of  ¥3  is 

svl =M+v 

8t  dt  2  dx2 


(0.28) 


(0.29) 


The  equations  for  ¥3  and  R12  should  be  the  same  because  the  quantities  are  equal  to  each  other  in  flows 
with  only  one  inhomogeneous  direction.  Therefore,  the  third  term  in  the  R)2  equation  must  appear  on  the 
right  hand  side  of  (0.27).  The  only  way  that  can  be  true  is  if  that  term  can  be  written  as  a  frame  invariant 

vector.  This  is  done  by  writing  R22Vj  2  as  2<|)S12 .  For  a  flow  with  only  one  direction  of  inhomogeneity  we 
replace  this  term  by  the  vector  -<j)s,  where  s  is  a  frame  invariant  term  to  be  defined  shortly.  In  this  way  we 
can  ensure  identical  frame  consistent  equations  for  R]2  and  in  flows  with  one  direction  of 
inhomogeneity,  where  they  are  equal  to  each  other. 
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